Using the basic ODE solver in Matlab
1. Example: Estimate the value of y(2), the solution to
y 00 + sin(t)y = 0,
y(−1) = 1, y 0 (−1) = 0
SOLUTION: We will actually do more and plot the solution as well. First we need an M-file for
the derivatives. Convert the second order D.E. into a system of first order, by defining u(t) = y(t),
v(t) = y 0 (t). Then:
u0 = v
v 0 = − sin(t)u
Now I’ll need a Matlab file for the derivatives:
function dy=Sample01(t,y)
% dy=Sample01(t,y) - A Sample DE for our ODE solver
dy=zeros(size(y));
dy(1)=y(2);
dy(2)=-sin(t)*y(1);
In Matlab, we type the commands:
tspan=linspace(-1,2);
%If we don’t want all this output, just do: tspan=[-1 2];
[tout,yout]=ode45(@Sample01,tspan,[1;0]);
plot(tout,yout);
The desire result is yout(end, 1) ≈ 0.6954
2. Example: Solve the Van der Pol Equation:
u00 − (1 − u2 ) + u = 0
with various initial conditions. Plot the resulting solution in the parameterized form (u(t), u0 (t)).
SOLUTION: In Matlab, we will need to write the differential equation as a first order D.E. Matlab
specifies that the function input (t, y) and output y 0 . We do this by introducing a “dummy” variable:
Let x = u, y = u0 . This gives a system of first order D.E.’s:
x0
y0
=y
= −x + (1 − x2 )y
I’ll write the following function and save it in my directory as vdp01.m:
function dy=vdp01(t,y)
%Function: dy=vdp(t,y) is the Van der Pol equation.
% Although t is never used, we still include it for the
% ODE solver.
dy=zeros(size(y));
dy(1)=y(2);
dy(2)=-y(1)+(1-y(1)^2)*y(2);
I’ll write the following as a script, but I could just type the commands in the command window:
1 tspan=linspace(0,40,200);
%I’ll have Matlab return the value
%my solution at these time values.
%I’ll use the following four initial conditions
%(could set these randomly)
IC=[2.5, 2.5; 2.5, -2.5; -2.1, 0.5; 0.5, 0.3];
%Now solve the DE four times
for j=1:4
[Tout,Yout{j}]=ode45(@vdp01,tspan,IC(j,:)’);
end
%Now plot the solutions:
figure(1)
for j=1:4
plot(Yout{j}(:,1),Yout{j}(:,1));
if j==1; hold on; end;
end
hold off
3. Example: Generate a solution to the Lorenz Equations (these are the “famous” values of the parameters):
x0 = 10(y − x)
y 0 = 28x − y − xz
z 0 = −(8/3) z + xy
I’ll once again write a function that takes in (t, y) and outputs the derivative:
function dy=lorenz01(t,y)
% dy=lorenz01(t,y) are the Lorenz Equations
dy=zeros(size(y));
dy(1)=10*(y(2)-y(1));
dy(2)=28*y(1)-y(2)-y(1)*y(3);
dy(3)=-(8/3)*y(3)+y(1)*y(2);
And I’ll use a script driver:
%% Script for the Lorenz Equations
tspan=linspace(0,75,4000); %Get lots of data!
IC=[5;5;5];
[Tout,Yout]=ode45(@lorenz01,tspan,IC);
plot3(Yout(:,1),Y(:,2),Y(:,3));
2