	Department of Engineering

University of Cambridge > Engineering Department > computing help > programs > matlab

Matlab and ODEs

With Matlab it's easy to solve some ODEs, though for harder ODEs you may need to do some preliminary maths first. This short document will illustrate the easier options, beginning with numerical solutions, then showing how to get general (symbolic) results.

A simple example

Suppose you want to solve

y' + 2ty =0
between t=0 and t=5
where y(0)=3

To be able to use the matlab routines you first need to write a function that calculates y' given t and y. That's easy to do here, which is why this example is easy.

function dy = yprime1(t,y)
dy= -2*t*y;

This function can go in a separate file called yprime1.m, though it's often tidier to make it a subfunction in the same file as the rest of the code (which is what we'll do here). Then all you need to do is run one of matlab's ODE-solver routines. If you put the following in a file called odedemo1.m and run it you'll see the result displayed in a graphics window.

function odedemo1
disp('Solving y''+2ty=0; y(0)=3; between t=0 and t=5');
[T,Y] = ode23(@yprime1, [0, 5], 3);
plot(T,Y);

function dy = yprime1(t,y)
dy= -2*t*y;

It's worth drawing the graph as a quick check even if you're only interested in the numbers. Note that ode23 in this example is given 3 arguments. The first defines the function (the @ before the name creates a "function handle" which is the recommended way of passing a function name). The 2nd and 3rd arguments describe the range and initial conditions respectively.

More arguments

Before trying harder ODEs, let's consider some variations on this simple example

More ODE options - To control the behaviour of the solver you can use odeset. For example, if you change odedemo1.m to

function odedemo1
disp('Solving y''+2ty=0; y(0)=3; between t=0 and t=5');
options=odeset('RelTol',1e-2);
[T,Y] = ode23(@yprime1, [0, 5], 3, options);
plot(T,Y);

function dy = yprime1(t,y)
dy= -2*t*y;

you'll get a relative error tolerance of 1% instead of the default 0.1%

Passing parameters to the function - You may wish to make the subfunction more flexible by letting it accept more arguments. For example, you may want to vary the constant in dy= -2*t*y. The following version of odedemo1.m compares results when using various parameter values, creating a new figure for each value.

function odedemo1
disp('Solving y''+pty=0; y(0)=3; between t=0 and t=5 for p=1, 2, and 3');

for i=1:3
   [T,Y] = ode23(@yprime1, [0, 5], 3, [], i);
   figure
   plot(T,Y);
end

function dy = yprime1(t,y, p)
dy= -p*t*y;

ode23

Second Order ODEs

In the following example we have a 2nd derivative

y'' + 2ty' +3y=0
between t=0 and t=5
where y(0)=3 and y'(0)=1

The trick is to reduce the order of this to produce 2 coupled 1st order ODEs by introducing some new variables

y₁ = y
y₂ = y'

Then the original system can be rewritten as

y₂' = -2ty₂ -3y₁
y₁'=y₂

ode23 can solve coupled systems like this, though vectors need to be used where in the earlier example we used single values. The subfunction (let's call it yprime2) now becomes

function dy = yprime2(t,y)
dy= [y(2); -2*t*y(2)-3*y(1)];

Note that matlab's y(1) and y(2) are what we called y₁ and y₂. The returned value is now a vector containing [y₁'; y₂'].

The initial conditions also need to be provided in a vector form. At t=0, [y₁; y₂] = [3;1]. Putting all this together we can create odedemo2.m

function odedemo2
[T,Y] = ode23(@yprime2, [0, 5], [3;1]);
plot(T,Y);

function dy = yprime2(t,y)
dy= [y(2); -2*t*y(2)-3*y(1)];

Now that ode23 is being given vectorised inputs it produces a vectorised output, so on the graph 2 lines will be displayed, showing y₁ and y₂ (which in terms of the original system are y and y').

More generally

the function provided to ode23 must return something that's the same size as the input parameter y
ode23 returns a vector of time values (in this case T - note that the values might not be uniformly spaced) and an array of solutions (in this case Y). Each column of Y corresponds to a different variable. There are as many rows as there are values in T.

Other ODE routines

Only ode23 has been used in these examples. It's not the only routine available, and it's not always the most accurate. If you have the time you could explore alternatives -

ode23 (non-stiff differential equations, low order method)
ode23s (stiff differential equations, low order method)
ode23t (non-stiff differential equations, low order method)
ode23tb (moderately stiff ODEs and DAEs, trapezoidal rule)
ode113 (non-stiff differential equations, variable order method)
ode15i (fully implicit differential equations, variable order method)
ode15s (stiff ODEs and DAEs Index 1, variable order method)
ode45 (non-stiff differential equations, medium order method)

General Solutions - the Symbolic toolbox

A matlab add-on called the Symbolic toolbox provides a routine called dsolve which gives general solutions. For example,

dsolve('D2y + 2*Dy - 3*y=0')

(where D is the differential operator) returns C1*exp(t)+C2*exp(-3*t) (where C1 and C2 are arbitrary constants). If boundary conditions are supplied, dsolve can give numerical solutions. E.g.

dsolve('D2y -5*Dy + 4*y=0', 'y(0) = 1', 'Dy(0)=0')

returns 4/3*exp(t)-1/3*exp(4*t)