Matlab/Octave

 

 

 

Differential Equation

 

 

ODE45

 

Lsim-State Space Model

 

Step

 

 

 

ode45 - 1st order

 

Ex)

Input

f = inline('-2*y','t','y');

odeopt = odeset ('RelTol', 0.00001, 'AbsTol', 0.00001,'InitialStep',0.5,'MaxStep',0.5);

[t,y] = ode45(f,[0 5],2,odeopt);

plot(t,y)

Output

 

 

 

ode45 - Undamped Pendulum

 

Ex 1)

Input

%Save the following contents in a .m file and run the .m file

omega = 1;

dy_dt = @(t,y) [y(2);...

                      -omega.^2*sin(y(1))];

odeopt = odeset ('RelTol', 0.00001, 'AbsTol', 0.00001,'InitialStep',0.5,'MaxStep',0.5);

[t,y] = ode45(dy_dt,[0 25], [0.0 1.0],odeopt);

subplot(1,2,1);plot(t,y(:,1),'r-',t,y(:,2),'b-');xlabel('time'); legend('y1(t)','y2(t)');

subplot(1,2,2);plot(y(:,1),y(:,2),'b-'); axis([-2.5 2.5 -2.5 2.5]); xlabel('y(2)');ylabel('y(1)');

Output

 

Ex 2) If you change the initial condition, you would get different result (a little bit distorted circle)

Input

%Save the following contents in a .m file and run the .m file

omega = 1;

dy_dt = @(t,y) [y(2);...

                      -omega.^2*sin(y(1))];

odeopt = odeset ('RelTol', 0.00001, 'AbsTol', 0.00001,'InitialStep',0.5,'MaxStep',0.5);

[t,y] = ode45(dy_dt,[0 25], [0.0 1.7],odeopt);

subplot(1,2,1);plot(t,y(:,1),'r-',t,y(:,2),'b-'); xlabel('time'); legend('y1(t)','y2(t)');

subplot(1,2,2);plot(y(:,1),y(:,2),'b-'); xlabel('y(2)');ylabel('y(1)');axis([-2.5 2.5 -2.5 2.5]);

Output

 

 

 

ode45 - Single Spring Mass- Undamped

 

 

 

Ex)

Input

%Save the following contents in a .m file and run the .m file

m = 1;

k = 1;

 

dy_dt = @(t,y) [y(2);...

                      -(k/m) * y(1) ];

odeopt = odeset ('RelTol', 0.00001, 'AbsTol', 0.00001,'InitialStep',0.5,'MaxStep',0.5);

[t,y] = ode45(dy_dt,[0 25], [0.0 1.0],odeopt);

subplot(1,2,1);plot(t,y(:,1),'r-',t,y(:,2),'b-'); xlabel('time'); ylim([-1.2 1.2]); legend('y1(t)','y2(t)');

subplot(1,2,2);plot(y(:,1),y(:,2),'b-'); xlabel('y(2)');ylabel('y(1)'); xlim([-1.2 1.2]); ylim([-1.2 1.2]);

Output

 

 

 

ode45 - Single Spring Mass- Damped

 

 

 

Ex)

Input

%Save the following contents in a .m file and run the .m file

m = 1;

k = 1;

c = 0.3;

 

dy_dt = @(t,y) [y(2);...

                 -(c/m) * y(2) - (k/m) * y(1) ];

odeopt = odeset ('RelTol', 0.00001, 'AbsTol', 0.00001,'InitialStep',0.5,'MaxStep',0.5);

[t,y] = ode45(dy_dt,[0 25], [0.0 1.0],odeopt);

subplot(1,2,1);plot(t,y(:,1),'r-',t,y(:,2),'b-'); xlabel('time'); ylim([-1.2 1.2]); legend('y1(t)','y2(t)');

subplot(1,2,2);plot(y(:,1),y(:,2),'b-'); xlabel('y(2)');ylabel('y(1)'); xlim([-1.2 1.2]); ylim([-1.2 1.2]);

Output

 

 

 

ode45 - Single Spring Mass- Damped and External Force

 

 

Ex)

Input

%Save the following contents in a .m file and run the .m file

m = 1;

k = 1;

c = 0.3;

F0 = 0.5;

w = 2.5;

 

dy_dt = @(t,y) [y(2);...

                      -(c/m) * y(2) - (k/m) * y(1) + (F0/m) * cos(w*t)];

odeopt = odeset ('RelTol', 0.00001, 'AbsTol', 0.00001,'InitialStep',0.5,'MaxStep',0.5);

[t,y] = ode45(dy_dt,[0 50], [0.0 1.0],odeopt);

subplot(1,2,1);plot(t,y(:,1),'r-',t,y(:,2),'b-'); xlabel('time'); ylim([-1.2 1.2]); legend('y1(t)','y2(t)');

subplot(1,2,2);plot(y(:,1),y(:,2),'b-'); xlabel('y(2)');ylabel('y(1)'); xlim([-1.2 1.2]); ylim([-1.2 1.2]);

Output

 

 

 

ode45 - Single Spring Mass- Damped and External Force with Frequency Sweep

 

Ex)

Input

%Save the following contents in a .m file and run the .m file

m = 1.0;

k = 0.7;

c = 0.3;

F0 = 0.5;

 

tmax = 100;

yinit = 0.0;

 

for i = 1:10

 

w = 0.2 * i;

dy_dt = @(t,y) [y(2);...

                 -(c/m) * y(2) - (k/m) * y(1) + (F0/m) * cos(w*t)];

odeopt = odeset ('RelTol', 0.00001, 'AbsTol', 0.00001,'InitialStep',0.5,'MaxStep',0.5);

[t,y] = ode45(dy_dt,[0 tmax], [0.0 yinit],odeopt);

subplot(10,7,[(i*7-6) (i*7-1)] );plot(t,y(:,1),'r-',t,y(:,2),'b-');

                                           ylim([-2 2]); set(gca,'xticklabel',[]);set(gca,'yticklabel',[]);

subplot(10,7,i*7);plot(y(:,1),y(:,2),'b-');

                        xlim([-2 2]); ylim([-2 2]); set(gca,'xticklabel',[]);set(gca,'yticklabel',[]);

 

end;

Output

 

 

 

ode45 - 1s Order System Equation- Lorenz Attractor

 

 

Ex)

Input

%Save the following contents in a .m file and run the .m file

P = 10;

r = 28;

b = 8/3;

 

dy_dt = @(t,y) [-P*y(1)+P*y(2);...

                    r.*y(1)-y(2)-y(1)*y(3);...

                    y(1)*y(2)-b.*y(3)];

odeopt = odeset ('RelTol', 0.00001, 'AbsTol', 0.00001,'InitialStep',0.5,'MaxStep',0.5);

[t,y] = ode45(dy_dt,[0 250], [1.0 1.0 1.0],odeopt);

subplot(1,3,1);plot(y(:,1),y(:,2),'r-'); xlabel('y(1)'); ylabel('y(2)');

subplot(1,3,2);plot(y(:,2),y(:,3),'g-'); xlabel('y(2)'); ylabel('y(3)');

subplot(1,3,3);plot(y(:,1),y(:,3),'b-'); xlabel('y(1)'); ylabel('y(3)');

Output

 

 

 

ode45 - Chemical Reaction

 

 

 

Ex)

Input

%Save the following contents in a .m file and run the .m file

k1 = 5;

k2 = 2;

k3 = 1;

 

dy_dt = @(t,y) [-k1*y(1) + k2*y(2);

                k1*y(1) - (k2+k3)*y(2);

                k3*y(2)];

odeopt = odeset ('RelTol', 0.00001, 'AbsTol', 0.00001,'InitialStep',0.5,'MaxStep',0.5);

[t,y] = ode45(dy_dt,[0 4], [1.0 0.0 0.0],odeopt);

plot(t,y(:,1),'r-',t,y(:,2),'g-',t,y(:,3),'b-'); xlabel('time'); legend('y(1)','y(2)','y(3)');

Output

 

 

 

Ode45 - Vander Pol Oscillator

 

 

 

 

 

Ex)

Input

%Save the following contents in a .m file and run the .m file

mu = 1;

k = 1;

 

dy_dt = @(t,y) [y(2);...

                      mu*(1-y(1)^2)*y(2)-k*y(1)];

odeopt = odeset ('RelTol', 0.00001, 'AbsTol', 0.00001,'InitialStep',0.5,'MaxStep',0.5);

[t,y] = ode45(dy_dt,[0 40], [1.0 2.0],odeopt);

subplot(1,3,[1 2]);plot(t,y(:,1),'r-',t,y(:,2),'g-'); xlabel('time'); legend('y(1)','y(2)');

subplot(1,3,3);plot(y(:,1),y(:,2)); xlabel('y(1)'); ylabel('y(2)');

Output

 

 

 

Ode45 - Vander Pol Oscillator with External Force

 

 

 

 

 

Ex)

Input

%Save the following contents in a .m file and run the .m file

mu = 1;

k = 1;

A = 2.0;

w = 10;

 

dy_dt = @(t,y) [y(2);...

                mu*(1-y(1)^2)*y(2)-k*y(1) + A*sin(w*t)];

odeopt = odeset ('RelTol', 0.00001, 'AbsTol', 0.00001,'InitialStep',0.5,'MaxStep',0.5);

[t,y] = ode45(dy_dt,[0 40], [1.0 2.0],odeopt);

subplot(1,3,[1 2]);plot(t,y(:,1),'r-',t,y(:,2),'g-'); xlabel('time'); legend('y(1)','y(2)');

subplot(1,3,3);plot(y(:,1),y(:,2)); xlabel('y(1)'); ylabel('y(2)');

Output

 

 

 

Ode45 - Duffing Oscillator

 

 

 

Ex)

Input

%Save the following contents in a .m file and run the .m file

delta = 0.06;

beta = 1.0;

w0 = 1.0;

w = 1.0;

gamma = 6.0;

phi = 0;

 

dy_dt = @(t,y) [y(2);...

                -delta*y(2)-(beta*y(1)^3 + w0^2*y(1))+gamma*cos(w*t+phi)];

odeopt = odeset ('RelTol', 0.00001, 'AbsTol', 0.00001,'InitialStep',0.5,'MaxStep',0.5);

[t,y] = ode45(dy_dt,[0 100], [3.0 4.1],odeopt);

subplot(1,3,[1 2]);plot(t,y(:,1),'r-',t,y(:,2),'g-'); xlabel('time'); legend('y(1)','y(2)');

subplot(1,3,3);plot(y(:,1),y(:,2)); xlabel('y(1)'); ylabel('y(2)');

Output

 

 

 

Ode45 - Matrix Equation - 2 x 2

 

 

 

 

Ex)

Input

%Save the following contents in a .m file and run the .m file

M = [-1,-1;...

     1,-2];

yinit = [1.0,1.0];

 

dy_dt = @(t,y) [M*y];

odeopt = odeset ('RelTol', 0.00001, 'AbsTol', 0.00001,'InitialStep',0.5,'MaxStep',0.5);

[t,y] = ode45(dy_dt,[0 4], yinit,odeopt);

plot(t,y(:,1),'r-',t,y(:,2),'g-'); xlabel('time'); legend('y(1)','y(2)');

Output

 

 

 

Lsim - Damped Spring

 

Ex) See Damped Spring Example in Differential Equation page for the description of the model.

 

Input

% This is tested on Octave. It would not work on Matlab (tested on Matlab 2014).

m = 1;

c = 0.2;

k = 1;

 

A = [0 1;-k/m -c/m];

B = [0 ; 1/m];

C = [0 1];

D = [0];

 

t = 0:0.1:50;

u = zeros(length(t),1);

 

x0 = [1 0];

 

sys = ss(A,B,C,0);

 

[y,t,x] = lsim(sys,u,t,x0);

plot(t,y);xlabel('time');ylabel('velocity');

Output

 

 

Ex)

 

Input

% This is tested on Octave. It would not work on Matlab (tested on Matlab 2014).

m = 1;

c = 0.2;

k = 1;

 

A = [0 1;-k/m -c/m];

B = [0 ; 1/m];

C = [1 0];

D = [0];

 

t = 0:0.1:50;

u = zeros(length(t),1);

 

x0 = [1 0];

 

sys = ss(A,B,C,0);

 

[y,t,x] = lsim(sys,u,t,x0);

plot(t,y);xlabel('time');ylabel('velocity');

Output

 

 

 

Lsim - RLC Circuit

 

Ex)  See RLC Circuit Example in Differential Equation page for the description of the model.

 

 

Input

% This is tested on Octave. It would not work on Matlab (tested on Matlab 2014).

R = 1;

L = 0.1;

C = 0.01;

 

A = [0 1;-1/(L*C) -R/L];

B = [0 ; 1/(L*C)];

C = [1 0];

D = [0];

 

t = 0:0.01:10;

u = ones(length(t),1);

u(200:400) = 0;

 

x0 = [1 0];

 

sys = ss(A,B,C,D);

 

[y,t,x] = lsim(sys,u,t,x0);

plot(t,y,'r-',t,u,'b--');xlabel('time');ylabel('voltage');legend('v2(t)','u(t)');

Output

 

 

 

Step - 1/s

 

 

Ex)

Input

sys = tf([1],[1 0]) % Create a transfer function based on numerator and denominator

step(sys);

Output

Transfer function 'sys' from input 'u1' to output ...

 

      1

 y1:  -

      s

 

Continuous-time model.

 

 

 

Step - 1/(a*s+1)

 

Ex)

Input

tau = 0.5;

sys = tf([1],[tau 1])

step(sys,5.0);

Output

Transfer function 'sys' from input 'u1' to output ...

 

          1

 y1:  ---------

      0.5 s + 1

 

Continuous-time model.

 

 

 

Step - 1/(a*s^2+b*s+1)

 

Ex)

Input

a = 1.5;

b = 0.7;

sys = tf([1],[a b 1])

step(sys,20.0);

Output

Transfer function 'sys' from input 'u1' to output ...

 

               1

 y1:  -------------------

      1.5 s^2 + 0.7 s + 1

 

Continuous-time model.