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:
- Group name.
- Name of the students in each group (with QUID #).
- Name of the project for each group.
You may check the selected projects by different groups using the link above.
• Written Report (Paper). Each group should submit 5-8 pages final report by
December 02, 2018. The report should be typed, should be of a professional quality
and should contain the following - Title-page.
- Abstract.
- Introduction.
- Methods.
- Results and Discussions.
- Conclusion.
- References.
- Any source codes and input/output files should be included as Appendix
• Authors’ Contributions: Each group should include a section at the end of the
paper called ”Authors’ Contributions”. In this section, you should report the contribution of each member in your group. If all members of the team fulfill their responsibilities correctly then each
member will receive the same points for the project. If
1
Math 217 Section L52
an individual’s portion of the work is incomplete or poorly completed then that individual will receive a different mark for the project. You may refer to any published
articles to know how this section can be written.
2 Specific Guidelines
The paper should include:
• A description of the mathematical problem which clearly identifies all the variables,
their type and dimensions, and any default values.
• A description of the differential equation model and the assumptions being made.
• A detailed discussion of the mathematical and/or numerical method used to solve the
problem. Include the derivations of formulas and use flow-charts and/or pseudocodes
to describe any numerical algorithms that may have been used.
• Appropriate discussion of the validation process used and an investigation of the
quality of the model.
3 Suggested Projects:
The following is a list of possible research projects from the textbook. These projects
focused on Engineering related problems. Chapter and page numbers refer to the course
textbook.
• Chapter 2 (pages 79-88):
- Oil Spill in a Canal
- Differential Equations in Clinical Medicine
- Torricellis Law of Fluid Flow
- The Snowplow Problem, Two Snowplows
- Clairaut Equations and Singular Solutions
- Multiple Solutions of a First-Order Initial Value Problem
- Utility Functions and Risk Aversion
- Designing a Solar Collector
- Asymptotic Behavior of Solutions to Linear Equations.
• Chapter 3 (pages 141-152): - Dynamics of HIV Infection
2
Math 217 Section L52 - Aquaculture
- Curve of Pursuit
- Aircraft Guidance in a Crosswind
- Bang-Bang Controls
- Market Equilibrium: Stability and Time Paths
- Period Doubling and Chaos
• Chapter 4 (pages 235-241): - Nonlinear Equations Solvable by First-Order Techniques
- Apollo Reentry
- Simple Pendulum
- Linearization of Nonlinear Problems
- Convolution Method
- Undetermined Coefficients Using Complex Arithmetic
- Asymptotic Behavior of Solutions
• Chapter 5 (pages 309-317): - Spread of Staph Infections in HospitalsPart I
- Hamiltonian Systems
- Cleaning Up the Great Lakes
- A Growth Model for PhytoplanktonPart I
• Chapter 8 (pages 495-497). - Buckling of a Tower
- Aging Spring and Bessel Functions
4 References
[1] Fundamentals of Differential Equations and Boundary Value Problems, by R.K. Nagle,
E.B. Saff and A.D. Snider, 8th Edition, 2012, Addison Wesly.
3
Solving ODEs
in Matlab
Dr. Houssein Ayoub
Department of Mathemtics, Statistics, and Physics, Qatar
University
- Outline –
I. Defining an ODE function in an M-file
II. Solving first-order ODEs
III. Solving systems of first-order ODEs
IV. Solving higher order ODEs
What are we doing when
numerically solving ODE’s?
Numerical methods are used to solve initial value
problems where it is difficult to obtain exact solutions
• An ODE is an equation that contains one independent variable (e.g. time)
and one or more derivatives with respect to that independent variable.
• In the time domain, ODEs are initial-value problems, so all the conditions
are specified at the initial time t = 0.
• Matlab has several different functions (built-ins) for the numerical
solution of ODEs. These solvers can be used with the following syntax:
[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
dt
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
!
*
*
*
*
- *
*
We know t0 and y(t0)
and we know the slope of
y(t), dy/dt =f(t,y).
We don’t know y(t) for any
values of t other than t0.
Use if ode45
fails because the
problem is stiff*
Low to
medium
ode15s
For
computationally
intensive
problems
ode113 Low to high
Less accurate
than ode45
ode23 Low
This should be the
first solver you try
ode45 Medium
Solver Accuracy Description
Runge-Kutta
(4,5) formula
*No precise definition of stiffness, but the main idea is that the equation
includes some terms that can lead to rapid variation in the solution.
[t,state]
= ode45(@dstate,tspan,ICs,options)
Defining an ODE function in an M-file
- Define tspan, ICs and options in one file (e.g.
call_dstate.m), which sets up ode45 - Define your constants and derivatives in another file
(e.g. dstate.m) or as a function dstate within the call
file - Run call_dstate.m
- Analyze the results as a plot of state vs. t
II. Solving first-order ODEs
!
dy
dt
= y'(t) =”y(t) # $y(t)
2
!
y(0) =10
Example:
function [t,y] = call_dstate()
tspan = [0 9]; % set time interval
y0 = 10; % set initial condition
% dstate evaluates r.h.s. of the ode
[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-