|
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.
|
|
|