Beginning Programming for
Numerical Solutions to Differential
Initial Value Problems
We wish to "solve" the ordinary differential equation (ODE):
where we have the initial conditions:
By "solve", we want to compute
Euler's Method
Matlab's ODE Solvers
Matlab provides a selection of solvers to integrate an ODE
function over a range. General form:
[t,y] = solver(odefun, tspan, y0)
The odefun is the function to be integrated over the interval
tspan with initial conditions y0. The result is the vector t for
values of the independent variable in the range tspan, and
estimated values y for each value t.
The solver is typically ode45.
Simple usage of ode45
We would like to solve this system numerically:
We need a function to calculate
function [dy] = simple_ode(t, y)
dy = -3.7*y;
Call the solver:
[t,y] = ode45(@simple_ode, [0 1], 5);
Using ode45, continued
Given the example:
[t,y] = ode45(@simple_ode, [0 1], 5);
• The interval [0 1] is the range over which the independent
variable (often x or t) varies. The solver chooses its own
values to use. Actual values chosen are sorted in t, a
column vector.
• Corresponding values of the integrated function are stored in
the column vector y.
• Instead of [0 1], additional points can be given, e.g,
0:0.1:1 -- which will be used as the values of the output t
(and to compute y).
Function Handles
The notation @simple_ode is a function handle, to allow the
solver to call our function simple_ode.
Small functions can be defined anonymously:
odefun = @(T,Y) -3.7*Y;
[t,y] = ode45(odefun, [0 1], 5);
or even:
[t,y] = ode45(@(T,Y)-3.7*Y, [0 1], 5);
Multiple ODEs
The Matlab ODE solvers can solve systems of related
equations over the same interval of the independent variable.
Example: The Lotka-Volterra predator-prey model can be used
to model the population on rabbits (r) and foxes (f):
Model this with
Rabbits and Foxes (2) ...
Need the ODE function:
function [dy] = foxrabbitode(t,y)
% FOXRABBITODE calculates population changes
% y(1) = rabbits
% y(2) = foxes
alpha = 0.01;
r = y(1);
f = y(2);
dy = [ 2*r - alpha*r*f
-f + alpha*r*f ];
• Calculate on y, a column vector, and produce a column
vector. Each row corresponds to one of the equations.
Rabbits and Foxes (3) ...
Create column vector of initial conditions (r0 and f0)
y0 = [ 300; 150 ];
[t,y] = ode45(@foxrabbitode, [0 20], y0);
figure; hold on;
plot(t, y(:,1), 'b'); % rabbits
plot(t, y(:,2), 'r'); % foxes
Each column of dy corresponds to the points computed for one
of the equations.
Function handles and arguments
We might want to pass alpha as an argument to the ODE
function [dy] = foxrabbitode2(t,y,alpha)
r = y(1);
f = y(2);
dy = [ 2*r - alpha*r*f
-f + alpha*r*f ];
Need anonymous function to pass this to ODE function:
odefun = @(T,Y) foxrabbitode2(T,Y,alpha);
[t,y] = ode45(odefun, interval, y0);
Higher-order ODEs
Treat higher-order ODEs as systems of related ODEs.
The ODE function receives both y and y', and computes y' and
y''. (Simply copy the received y' value!)
We need initial conditions for y and y'.
Consider solving:
Let: y(1) = 0, y'(1) = 1
Solve for y(t), for t in [1 30].
Function for simple second-order ODE
function [dy] = second_ode(t,y)
% SECOND_ODE is a simple second order ODE.
% y(1) = f(t) % 1st equation
% y(2) = f'(t) % 2nd equation
% dy(1) = f'(t) % 1st equation
% dy(2) = f''(t) % 2nd equation
% Note how we can just copy y(2) to dy(1).
dy = [ y(2)
(-1/t)*sin(y(1)) ];
Solving the second-order ODE
% Compute the solutions.
[t,y] = ode45(@second_ode, [1 30], [0;1]);
% y(:,1) is y(t), and y(:,2) is y'(t).
plot(t, y(:,1))
Trajectory as ODE
function [dpv] = update_path(t, pv)
% pv = position and velocity state vector
pv(1) = x position, xpos, at time t
pv(2) = y position, ypos, at time t
pv(3) = x velocity, vx, at time t
pv(4) = y velocity, vy, at time t
x = pv(1);
y = pv(2);
xv = pv(3);
yv = pv(4);
dpv = [xv
% dx = x velocity
% dy = y velocity
% d2x = d(x velocity) = x acceleration
-9.81 ]; % d2y = d(y velocity) = y acceleration
Event functions
Sometimes we need to stop the integration due to some event,
especially if we don't know how long to run!
You can create an "event function" to stop integration:
function [value,isterminal,direction] = event(t,y)
This accepts t,y like the ODE function. It computes a value
at t,y. The output isterminal, if 1, means to stop the
integration if value is 0. The direction indicates when we
care about value being 0; if it is 0, then we always care.
Event functions: Trajectory
Sample for testing if above ground:
function [value, isterminal, direction] = ...
trajectory_event(t, pv)
if pv(2) > 0
value = 1; % Above the ground
value = 0; % Hit the ground
isterminal = 1; % Stop when value=0
direction = 0; % Always care about zeros
Using an event function
Use odeset to configure options to the ODE solvers, including
events functions:
options = odeset('Event', @trajectory_event);
[t,y] = ode45(@update_path, 0:dt:ft, ...
[x;y;v0x;v0y], options);
• You can use [0 Inf] as an interval when you don't know
when to finish, but have n Event function to stop.
• You can use more complicated function handles to pass
information to the Event function, e.g., the terrain.
Other solvers
• Some problems are stiff, needing slower specialized
solvers. Other solvers are faster than ode45, but not as
• Usually try ode45 first. If it works but is too slow, consider
• If ode45 is so slow it is unusable, you might have a stiff
problem. Try using ode15s or ode23s.
• See doc ode45 for much more information, including
guidance about choosing a solver, options you can pass to
solvers, etc.

similar documents