Math 217 Section L52
Project Guidelines
Dr. Houssein Ayoub
Department of Mathematics, Statistics, and Physics, Qatar University, Doha, Qatar
1 General Guidelines
The main objective of the project is to give you the opportunity to further develop the differential equation theories developed in the class and apply them to solve some interesting
differential equation models. A list of projects is suggested in the textbook (See Section
3). A project may be completed by a group of four to five students maximum. Please note
the following deadline and instructions.
• Team Selection and Project Selection. You should decide your teammates and
choose the project from the textbook by October 31, 2018. Please note that it
is not allowed to choose the same project by different groups. Each group
should fill the following information by the date above using the link here:
Solving ODEs
in Matlab
Dr. Houssein Ayoub
Department of Mathemtics, Statistics, and Physics, Qatar
University
[outputs]
= function_handle(inputs)
[t,state]
= solver(@dstate,tspan,ICs,options)
Matlab algorithm
(e.g., ode45,
ode23)
Handle for function
containing the
derivatives
Vector that specifiecs the
interval of the solution
(e.g., [t0:5:tf])
A vector of the
initial conditions
for the system
(row or column)
An array. The solution of
the ODE (the values of
the state at every time).
!
dy
t
y
!
y(0) =1
!
y(t) = t
2 +1
What are we doing when
numerically solving ODE’s?
Integrators compute nearby value of y(t) using
what we already know and repeat.
Higher order numerical methods reduce error at the cost of speed:
• Euler’s Method – 1st order expansion
• Midpoint method – 2nd order expansion
• Runge-Kutta – 4th order expansion
t t 0
y(t)
y(t0)
y
!
*
*
*
*
[t,state]
= ode45(@dstate,tspan,ICs,options)
Defining an ODE function in an M-file
[t,y]
= ode45( @dstate ,tspan ,y0);
plot(t,y)
disp([t,y]) % displays t and y(t)
function dydt = dstate (t,y)
alpha=2; gamma=0.0001;
dydt = alpha* y-gamma y^2;
end
end
1
2-
3-
4
5-
6-
7-
8
9-
10-
11-
12-
Save as call_dstate.m in some directory, and cd to that directory in the matlab GUI
Matlab ode45’s numerical solution
At t = 9, have we reached
!
steady state?
dy
dt
= y'(t) =”y(t) # $y(t)
2
!
y(0) =10
!
limt”># y(t) =
$
%
= 20,000 EDU>> [t, y] = call_dstate;
EDU>> steady_state = y(end)
steady_state =
1.9999e+04
From the command line:
III. Solving systems of first-order ODEs
• This is a system of ODEs because we have more than one derivative with
respect to our independent variable, time.
• This is a stiff system because the limit cycle has portions where the
solution components change slowly alternating with regions of very sharp
change – so we will need ode15s.
• This is a example from mathworks, a great resource @ mathworks.com or
the software manual.
• This time we’ll create separate files for the call function (call_osc.m) and
the ode function (osc.m)
!
dy1
dt
= y2
dy2
dt
=1000(1″ y1
2
)y2
” y1
!
y1(0) = 0
y2 (0) =1
van der Pol equations in relaxation oscillation:
III. Solving systems of first-order ODEs
!
dy1
dt
= y2
dy2
dt
=1000(1″ y1
2
)y2
” y1
!
y1(0) = 0
y2 (0) =1
van der Pol equations in relaxation oscillation:
To simulate this system, create a function osc containing the equations.
Method 1: preallocate space in a column vector, and fill with derivative functions
function dydt = osc(t,y)
dydt = zeros(2,1); % this creates an empty column
%vector that you can fill with your two derivatives:
dydt(1) = y(2);
dydt(2) = 1000(1 – y(1)^2)y(2) – y(1);
%In this case, y(1) is y1 and y(2) is y2, and dydt(1)
%is dy1/dt and dydt(2) is dy2/dt.
end
1
2-
3
4-
5-
6
7
8-
Save as osc.m in the same directory as before.
III. Solving systems of first-order ODEs
!
dy1
dt
= y2
dy2
dt
=1000(1″ y1
2
)y2
” y1
!
y1(0) = 0
y2 (0) =1
van der Pol equations in relaxation oscillation:
function dydt = osc(t,y)
dydt = [y(2)
1000(1 – y(1)^2)y(2) – y(1)];
%Still y(1) is y1 and y(2) is y2, and dydt(1)
%is dy1/dt and dydt(2) is dy2/dt.
end
1
2-
3
4
5
6-
Save as osc.m in the same directory as before.
To simulate this system, create a function osc containing the equations.
Method 2: vectorize the differential functions
III. Solving systems of first-order ODEs
!
dy1
dt
= y2
dy 2
dt
=1000(1″ y1
2
)y2
” y1
!
y1(0) = 2
y2(0) = 0
van der Pol equations in relaxation oscillation:
1
2-
3-
4-
5-
6-
7-
Save as call_osc.m in the same directory as before.
Now solve on a time interval from 0 to 3000 with the above initial conditions.
Create a scatter plot of y1 with time.
function [T,Y] = call_osc()
tspan = [0 3000];
y1_0 = 2;
y2_0 = 0;
[T,Y] = ode15s(@osc,tspan,[y1_0 y2_0]);
plot(T,Y(:,1),’o’)
end
!
dy1
dt
= y2
dy2
dt
=1000(1″ y1
2
)y2
” y1
!
y1(0) = 2
y2(0) = 0
van der Pol equations in
relaxation oscillation:
Plot of numerical solution
IV. Solving higher order ODEs
• Second order non-linear ODE
• Convert the 2nd order ODE to standard form:
2
2
2
2
sin
sin
d ML MG
dt
d G
dt L
! !
! !
= ”
= ”
MG
Simple pendulum:
!
z1 = “, z2 = d” /dt
dz1
dt
= z2
dz2
dt
= #
G
L sin(z1)
Non-linear pendulum function file
• G = 9.8 m/s
• L = 2 m
• Time 0 to 2π
• Initial θ = π/3
• Try ode23
• Plot θ with time
!
z1 = “, z2 = d” /dt
dz1
dt
= z2
dz2
dt
= #
G
L sin(z1)
function [] = call_pend()
tspan=[0 2pi]; % set time interval
z0=[pi/3,0]; % set initial conditions
[t,z]
=ode23(@pend,tspan,z0);
plot(t,z(:,1))
function dzdt = pend(t,z)
G=9.8; L=2; % set constants
z1=z(1); % get z1
z2=z(2); % get z2
dzdt = [z2 ; -G/L*sin(z1);];
end
end
1
2-
3-
4-
5-
6
7-
8-
9-
10-
11-
12-