Ordinary Differential Equations (ODE) in MATLAB Ordinary Differential Equations (ODE) in MATLAB

What will we learn from the next 5 lectures

Ordinary Differential Equations (ODE) in

MATLAB ◮ How to solve ODEs using MATLAB.
◮ How to model biological systems using ODEs in MATLAB.
Shan He ◮ How to analyse ODEs using MATLAB.
School for Computational Science
◮ Understand bifurcation and chaos using MATLAB.
University of Birmingham ◮ Applications of bifurcation and chaos to biological problems.

Module 06-23836: Computational Modelling with MATLAB

Ordinary Differential Equations (ODE) in MATLAB Ordinary Differential Equations (ODE) in MATLAB
Outline Concepts about ODE

Outline of Topics What is ODE

Concepts about ODE

An Ordinary Differential Equation (ODE) is an equation involving

Solving ODE in MATLAB
a function and its derivatives.
ODE Solvers in MATLAB
Solving linear ODEs in MATLAB
Solving high order ODEs in MATLAB
Solving ODEs in MATLAB: Advanced topics

Ordinary Differential Equations (ODE) in MATLAB Ordinary Differential Equations (ODE) in MATLAB
Concepts about ODE Concepts about ODE

Definition of ODE Linear ODE and Homogeneous Linear ODE

◮ A ODE is said to be linear if F can be written as a linear
combination of the derivatives of y together with a constant
◮ Let y be an unknown function of x. term, all possibly depending on x:
dy
◮ Let y ′ = dx be the first derivative with respect to x.
◮ Let y (n) = dny
be the nth derivative with respect to x. an (x)y n + an−1 (x)y n−1 + · · · + a1 (x)y ′ + a0 y = r (x)
dx n
◮ Let F be a given function or more concisely,
◮ Then an ODE of order n is an equation of the form: n−1
X
yn = ai (x)y (i) + r (x)
F (x, y , y ′ , . . . , y (n) ) = 0 i=0

where ai (x) and r (x) are continuous functions in x and the

function r (x) is called the source term.
◮ A linear ODE is said to be homogeneous if r (x) = 0,
otherwise it is called non-homogeneous or inhomogeneous.
Ordinary Differential Equations (ODE) in MATLAB Ordinary Differential Equations (ODE) in MATLAB
Solving ODE in MATLAB Solving ODE in MATLAB
ODE Solvers in MATLAB ODE Solvers in MATLAB

Solution to ODE ODE Solvers in MATLAB

◮ Matlab has several different ODE solvers for the numerical

◮ If an ODE is linear, it can be solved by analytical methods. solution of ODEs:
◮ In general, an nth-order ODE has n linearly independent ◮ ode45: based on an explicit Runge-Kutta (4, 5) formula and
solutions. the Dormand-Prince method.
◮ Any linear combination of linearly independent functions ◮ ode23: based on an explicit Runge-Kutta (2, 3) formula and
solutions is also a solution. the Bogacki and Shampine method.
◮ If an ODE is not linear, most of the case it cannot be solved ◮ We choose according to order of accuracy and the type of
exactly: we will use MATLAB to obtain approximate solutions systems (stiff or nonstiff).
◮ Rule of thumb: Always try ode45 first.

Ordinary Differential Equations (ODE) in MATLAB Ordinary Differential Equations (ODE) in MATLAB
Solving ODE in MATLAB Solving ODE in MATLAB
ODE Solvers in MATLAB Solving linear ODEs in MATLAB

How to use MATLAB ODE Solvers Solving linear ODE in MATLAB

◮ The MATLAB ODE solvers can be called as a function:

[T,Y] = ode**(@odefun,tspan,y0,options) One simple example:
◮ @odefun: Function handle of the ODE function (
du(x)
u′ = dx = u(x) + v (x)
◮ tspan: Interval of integration, [t0 ,tfinal ]. dv (x)
v′ = dx = u(x)
◮ y0: Initial conditions.
◮ options: Optional parameters. Analytical solution can be obtained, but how to solve it in
MATLAB numerically?
◮ The ODE function odefun define the ODEs:
function [dy] = odefun(t, y, parameters)

Ordinary Differential Equations (ODE) in MATLAB Ordinary Differential Equations (ODE) in MATLAB
Solving ODE in MATLAB Solving ODE in MATLAB
Solving linear ODEs in MATLAB Solving linear ODEs in MATLAB

Solving linear ODE in MATLAB Solving linear ODE in MATLAB

Steps for solving this equation numerically in MATLAB:

◮ Step 2: Call a numerical solver provided in MATLAB, e.g.,
◮ Step 1: Create a MATLAB function that defines the ODEs [T,Y] = ode45(odefun,tspan,y0)
1 function dy = simple_ode(t,y)
2 % function to be integrated 1 tspan = [0 10];
3 dy = zeros(2,1); 2 y0 = [1 0];
4 dy(1) = y(1) + y(2); 3 [X,Y] = ode45(@simple_ode,tspan,y0);
5 dy(2) = y(1);
Ordinary Differential Equations (ODE) in MATLAB Ordinary Differential Equations (ODE) in MATLAB
Solving ODE in MATLAB Solving ODE in MATLAB
Solving high order ODEs in MATLAB Solving high order ODEs in MATLAB

Reduction of ODE order Reduction of ODE order

Methods:
◮ Recall an ODE of the general form: ◮ We will use a second order ODE as an example:
′ (n)
F (x, y , y , . . . , y )=0
(
x ′ = −ye (−t/5) + y ′ e (−t/5) + 1
◮ We said this system is an ODE of order n. y ′′ = −2 sin(t)
◮ Any differential equation of order n can be reduced to a ◮ Step 1: Introduce a new variable that equals the first
system of n first-order (n = 1) differential equations. derivative of the free variable in the second order equation:
◮ We do so because high order ODE (n > 1) is difficult to solve. z = y′
◮ MATLAB ODE solvers cannot handle higher order ODE! ◮ Step 2: Taking the derivative of each side yields the following:
z ′ = y ′′

Ordinary Differential Equations (ODE) in MATLAB Ordinary Differential Equations (ODE) in MATLAB
Solving ODE in MATLAB Solving ODE in MATLAB
Solving high order ODEs in MATLAB Solving high order ODEs in MATLAB

Reduction of ODE order Solve high order ODE in MATLAB

Methods:
◮ Step 4: Write a MATLAB function ho ode.m to define the
Methods:
ODE:
◮ Step 3: Substituting the second order ODE with z and z ′ : 1 function dy = high_order_ode_example(t,x)

′ (−t/5) + ze (−t/5) + 1
2 % x(1) = x
x = −ye
 3 % x(2) = y
′
z = −2 sin(t)
 4 % x(3) = z
 ′ dy = [-x(2) * exp(-t/5) + ...
y =z 5
6 x(3) * exp(-t/5) + 1;
7 x(3);
8 -2*sin(t)]
9 end

Ordinary Differential Equations (ODE) in MATLAB Ordinary Differential Equations (ODE) in MATLAB
Solving ODE in MATLAB Solving ODE in MATLAB
Solving high order ODEs in MATLAB Solving ODEs in MATLAB: Advanced topics

Solve high order ODE in MATLAB Stiffness of ODE equations

◮ Stiffness is a subtle, difficult, and important concept in the
numerical solution of ordinary differential equations.
Methods: ◮ “A problem is stiff if the solution being sought varies slowly,
◮ Step 5: evaluate the system of equations using ODE45: but there are nearby solutions that vary rapidly, so the
1 t0 = 5; % Start time numerical method must take small steps to obtain satisfactory
2 tf = 20; % Stop time results.”
3 x0 = [1 -1 3] % Initial conditions ◮ Example:
4 [t,s] = ode45(@ho_ode,[t0,tf],x0);
5 x = s(:,1); y′ = y2 − y3
6 y = s(:,2);
7 z = s(:,3); y (0) = τ
8 plot(t,s); 2
0≤t≤
τ
◮ If we weren’t concerned with efficiency of ODE solver, we
wouldn’t be concerned about stiffness.
Ordinary Differential Equations (ODE) in MATLAB Ordinary Differential Equations (ODE) in MATLAB
Solving ODE in MATLAB Solving ODE in MATLAB
Solving ODEs in MATLAB: Advanced topics Solving ODEs in MATLAB: Advanced topics

Passing parameters to ODEs Passing parameters to ODEs

◮ Call the ODE solver and pass the parameter to the ODE
function:
◮ By passing the value to ODE function, we can change the 1 tspan = [0 10];
coefficients in ODEs each time we call the ODE solver. 2 initvalue = [1 0];
◮ We need to define the parameters in the ODE function header: 3 option = [];
1 function dy = ode_pv(t,y,A) 4 for ii=-10:4
2 % function to be integrated 5 A = [1 ii; 2 -1];
3 dy = zeros(2,1); 6 [t,y] = ode45(@ode_pv,...
4 dy = A*y 7 tspan,initvalue,option,A);
8 hold on;
9 plot(t, y);
10 end

Ordinary Differential Equations (ODE) in MATLAB Ordinary Differential Equations (ODE) in MATLAB
Solving ODE in MATLAB Solving ODE in MATLAB
Solving ODEs in MATLAB: Advanced topics Solving ODEs in MATLAB: Advanced topics

Events Define events

◮ Recall syntax of the ODE solvers:
[T,Y] = ode**(@odefun,tspan,y0,options) ◮ Define the ODE function:
◮ We generally assume tspan is known, e.g., 1 function dydt = ode_ball(t,y)
2 dydt = [y(2); -9.8];
t0 ≤ t ≤ tfinal ◮ Define event function for this problem:
1 function [value,isterminal,direction] ...
◮ But sometimes it is also important to determine tfinal .
2 = events(t,y)
◮ Example: a ball is falling because of gravity, when does it hit
3 % Locate the time when height passes through
the ground?
4 % zero and stop integration.
◮ Equations: 5 value = y(1); % Detect height = 0
6 isterminal = 1; % Stop the integration
y ′′ = g = −9.8 7 direction = -1; % Desired directionality of
where y (t) is the height of the object at time t 8 % the zero crossings:
◮ The question: for what t does y (t) = 0? We can use events.

Ordinary Differential Equations (ODE) in MATLAB Ordinary Differential Equations (ODE) in MATLAB
Solving ODE in MATLAB Solving ODE in MATLAB
Solving ODEs in MATLAB: Advanced topics Solving ODEs in MATLAB: Advanced topics

Run ODE Solver with events Other ODE solvers

◮ Define the event function as an option for ODE solver: ◮ What we have introduced are all for Initial Value Problems for
1 opts = odeset(’events’,@ball_events); ODEs.
◮ Solve it using MATLAB ode45: ◮ To solve Boundary Value Problems: bvp4c.
2 y0 = [20; 0]; Tutorial.
3 [t,y,tfinal] = ode45(@ode_ball,[0 Inf]... ◮ To solve Delay ODEs: dde23.
4 ,y0,opts); Tutorial.
5 tfinal ◮ To solve Stochastic ODEs: MATLAB SDE Toolbox.
6 plot(t,y(:,1),’-’,[0 tfinal],[1 0],’o’) Tutorial.