> For the complete documentation index, see [llms.txt](https://ykkim.gitbook.io/EC/llms.txt). Markdown versions of documentation pages are available by appending `.md` to page URLs; this page is available as [Markdown](https://ykkim.gitbook.io/EC/numerical-programming/matlab-example.md).

# Example: MATLAB

## MATLAB Online Course

Learn the basics of MATLAB® through this introductory tutorial on commonly used features and workflows.

{% embed url="<https://matlabacademy.mathworks.com/details/matlab-onramp/gettingstarted>" %}

List of other Self-Paced Online Courses:

{% embed url="<https://matlabacademy.mathworks.com/?page=1&sort=featured>" %}

## MATLAB Examples of Numerical Programming

by Y.-K. Kim

Example codes using MATLAB built-in functions

[You can download matlab code](https://github.com/ykkimhgu/NumericalProg-student/blob/main/tutorial/NP%20MATLAB%20Example%20codes.mlx)

### Non-linear solver

First, define a non-linear function in the form of f(v)=0.

Then, use a non-linear equation solver fzero (). ![image-20230818144953414](https://github.com/ykkimhgu/course-doc/assets/38373000/d6167efa-b2b9-4098-86be-68b2cc108cb9)

​

<figure><img src="/files/LC40TKoOpqhPHDmKGlYk" alt="" width="375"><figcaption></figcaption></figure>

```matlab
v0=0.5;
v=fzero(@fnSolarCell,v0)

function y=fnSolarCell(x)
    q=1.6022E-19; k=1.3806E-23; Voc=0.5; T=297;
    qkT=q/(k*T);
    y=exp(qkT*x)*(1+qkT*x)-exp(qkT*Voc);
end

```

### Integration

#### Integrating discrete dataset: trapz()

Trapezoidal Method

```matlab
% Discrete dataset
x=[0 5 10 15 20 25 30 35 40 45 50 55 60]
y=[0 3 8 20 33 42 40 48 60 12 8 4 3]
plot(x,y)

% Matlab function
I_matlab = trapz(x,y);     
```

#### Integrating a Function: integral(fun, a,b)

The area of the shaded region shown in the figure can be calculated by:

![image](https://github.com/ykkimhgu/course-doc/assets/38373000/4d64fc21-25cd-4346-aebc-10f8060612cd)

<figure><img src="/files/pOsxzniT3zWwj69V3ZnF" alt=""><figcaption></figcaption></figure>

```matlab
a=-3; b=3; N=8;
h=(b-a)/N;
x=a:h:b;
fun=@(x) (1-x.^2).^0.5;
% Matlab function
I_matlab = integral(fun,a,b);       
```

### Differentiation

#### Differentiating discrete dataset

```matlab
% Differentiation from discrete data
X = [1 1 2 3 5 8 13 21];
Y = diff(X)

% Differentiation from discrete data
h = 0.001;       % step size
X = -pi:h:pi;    % domain
f = sin(X);      % range
Y = diff(f)/h;   % first derivative
Z = diff(Y)/h;   % second derivative
plot(X(:,1:length(Y)),Y,'r',X,f,'b', X(:,1:length(Z)),Z,'k')
```

#### Differentiate a Function

```matlab
% Ordinary Differentiation of a function
syms x
g = exp(x)*cos(x);
diff(g)

% 2nd order Differentiation of a function
diff(g,2)

% Partial Differentiation of a function
syms s t
f = sin(s*t);
diff(f,t)
```

### Linear Equations

Solve for Ax=b

<figure><img src="/files/xpqXo5Ap5l405IXF3psZ" alt=""><figcaption></figcaption></figure>

![image](https://github.com/ykkimhgu/course-doc/assets/38373000/9e13e283-77a7-4454-9b32-58a80fa822fc)

```matlab
% A, b
A=[9.5, -2.5, 0, -2, 0;     -2.5, 11, -3.5, 0, -5;     0,-3.5, 15.5, 0, -4;     -2,  0,  0, 7, -3;     0, -5,  -4, -3, 12];
b=[12; -16; 14; 10; -30];

% solve for Ax=b
x=A\b   
x=inv(A)*b
% condition number
c=cond(A)
% norm 
n=norm(A)
% eigenvalue/vector
[eigVec,eigVa]=eig(A)
% QR factorization
[Q,R]=qr(A)

% LU factorization
[L,U]=lu(A)   
```

### Polynomial fitting

```matlab
t = 1:1:15;
V=[2.272	2.092	1.887	1.629	1.482	1.308	1.030	0.875	0.693	0.470	0.336	0.095	-0.163	-0.371	-0.511];

% Matlab function for polynomial fit
Zopt=polyfit(t,V,1);
Yopt=polyval(Zopt,t);  % Matlab function

figure
plot(t,V, '*r')
hold on
plot(t,Yopt, '-b')
xlabel('time','fontsize',15)
ylabel('V','fontsize',15)
```

### Exponential Fitting

RC circuit with unknown capacitor C and resistor of 5M

a) Find the capacitance C from curve fitting b) Estimate the voltage when time=32sec

<figure><img src="/files/7ARyeWdHMah9rAfANZxu" alt=""><figcaption></figcaption></figure>

![image](https://github.com/ykkimhgu/course-doc/assets/38373000/92ee0d7d-9180-4ddd-a2f9-58319afd0290)

```matlab
Xdata = 1:1:15;
Ydata = [9.7 8.1 6.6 5.1 4.4 3.7 2.8 2.4 2.0 1.6 1.4 1.1 0.85 0.69 0.6];


% Matlab function
Zopt=polyfit(Xdata,log(Ydata),1)
R=5e6;
a0=Zopt(2);
a1=Zopt(1);
V0=exp(a0);
tau=-1/a1;
C=tau/R;

% Exponential model
Yopt=V0*exp(-1/(R*C).*Xdata);

figure
plot(Xdata,Ydata, '*r')
hold on
plot(Xdata,Yopt)
hold on
xlabel('time (s)','fontsize',15)
ylabel('V_R(volt)','fontsize',15)
```

### 1st order ODE-IVP

Solve for the response of an RC circuit with a DC input. Let tau= 4.9919; Vm=11.91;

<figure><img src="/files/nHN6us97r6pMCB9g9hnE" alt=""><figcaption></figcaption></figure>

![image](https://github.com/ykkimhgu/course-doc/assets/38373000/cc43d271-fdd5-4d27-ab8f-f0959338d81c)

```matlab
% Initial Condition
a=0; b=15; h=1; 
y0 = 11.91;
t=a:h:b;

%% MATLAB's function ODE45
[tmat,vmat] = ode45(@myRC, [a b], y0);  % Fourth/Fifth RK
figure()
plot(tmat,vmat,'.b')
xlabel('time (s)','fontsize',15)
ylabel('V_R(volt)','fontsize',15)

function dvdx = myRC(v)
    tau=4.9919; T=1/tau;  Vm=11.91;
    dvdx =-T*v + 1*T*Vm;
end

```

### 2nd order ODE - IVP

Solve an m-c-k system with a sinusoidal input.

Use m=10kg ; k=800 N/m; c=200 N/(m/s), f=10Hz , h=0.01, Fdc=100N.

<figure><img src="/files/1GD0u7NCDqkYUR6h45Q7" alt=""><figcaption></figcaption></figure>

![image](https://github.com/ykkimhgu/course-doc/assets/38373000/353efdc1-7ba3-4a47-9931-cee8617b067d)

```matlab
% Initial Condition
y0 = 0; v0 = 0;
Yinit = [y0 v0];
a=0; b=1; h=0.01;
tspan = [a:h:b];

% MATLAB's function ODE45
[Time Y] = ode45(@mckFunc,tspan,Yinit);

figure
subplot(2,1,1)
    plot(Time,Y(:,1),'--b')
    xlabel('Time (s)')
    ylabel('Position (m)')
    title('ode45')
subplot(2,1,2)
    plot(Time,Y(:,2),'k')
    xlabel('Time (s)')
    ylabel('Velocity (m/s)')

function [dXdt] = mckFunc(t,x)
    dXdt=zeros(2,1); % column vector
    m=10; k=800; c=200; f=10;
    FinDC=100;
    Fin=FinDC*cos(2*pi*f*t);
    dXdt(1)=x(2);
    dXdt(2)=1/m*(Fin-c*x(2)-k*x(1));
end
```

### Eigenvalue

What are the eigenvalues for a given m-c-k system response? Use m=10kg ; k=800 N/m; c=200 N/(m/s)

<figure><img src="/files/PC1rpSarFvRMAgqEJCyg" alt="" width="174"><figcaption></figcaption></figure>

<figure><img src="/files/hmilQAuk9r0pkcNWtfBV" alt="" width="375"><figcaption></figcaption></figure>

```matlab
k=800; c=200; m=10;
A = [0 1; -k/m, -c/m];
disp('Eigvalue and vector of A (MATLAB):');
[eigVec,eigVa]=eig(A)

```

## More tutorial codes

### [NP lecture tutorial codes](https://github.com/ykkimhgu/NumericalProg-student/tree/main/tutorial)


---

# Agent Instructions
This documentation is published with GitBook. GitBook is the documentation platform designed so that both humans and AI agents can read, navigate, and reason over technical content effectively. Learn more at gitbook.com.

## Querying This Documentation
If you need additional information that is not directly available in this page, you can query the documentation dynamically by asking a question.

Perform an HTTP GET request on the current page URL with the `ask` query parameter:

```
GET https://ykkim.gitbook.io/EC/numerical-programming/matlab-example.md?ask=<question>
```

The question should be specific, self-contained, and written in natural language.
The response will contain a direct answer to the question and relevant excerpts and sources from the documentation.

Use this mechanism when the answer is not explicitly present in the current page, you need clarification or additional context, or you want to retrieve related documentation sections.