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
you'll get a relative error tolerance of 1% instead of the default 0.1%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; - 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;
The empty square brackets in the call to ode23 are a place-holder because this
time we don't want to add any options.
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,/p>
y_{1} = y y_{2} = y' |
Then the original system can be rewritten as
y_{2}' = -2ty_{2} -3y_{1} y_{1}'=y_{2} |
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_{1} and y_{2}. The returned value is now a vector containing [y_{1}'; y_{2}'].
The initial conditions also need to be provided in a vector form. At t=0, [y_{1}; y_{2}] = [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_{1} and y_{2} (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)
See Also
Whole books exist on the topic of using matlab to solve ODEs. Here are just a few of the online resources.
- Matlab ODE Suite (a Mathworks Technical Paper)
- Using Matlab for First Order ODEs (Tobias von Petersdorff)
- Using Matlab for Higher Order ODEs and Systems of ODEs (Tobias von Petersdorff)
- Solving ODEs with Matlab (M.A. Cappelli). This document provided some of the examples used on this page
- Using ODE45 (Bucknell University)