MATLAB Programming/Advanced Topics/Numerical Manipulation/Ordinary Differential Equations

All of the differential equations have the same syntax that you must use, and the same input and output arguments. All of the solvers require three input arguments: a function handle to the differential equation you want to solve, a time interval over which to integrate, and a set of initial conditions. Let us suppose we want to solve the simple differential equation y' = y, y(0) = 1, which has the true solution y(t) = e^t. Suppose we want the solution between t = 0 and t = 1.

To use function handles, you must first create an M-file with the function in it like so:

function Yprime = func(t, y)
Yprime = 0;
Yprime = y;

Note that you must include the time argument even if it is not used in the differential equation. The initialization of Yprime to an array of zeros will save you grief if you try to solve more than one function; Yprime must be returned as a VERTICAL array but, if you don't initialize it as a verical array (or transpose at the end), it will return a HORIZONTAL array by default.

Also note that, with the exception of ode15i, the function must be solved explicitly for y'.

Once this file is created, call the ODE function with the arguments in the order (function, timeinterval, initialcond):

>> ode45(@func, [0,1], 1)

The other method is to use anonymous functions, which is only useful if you have one function (otherwise you must use function handles). You again must declare the anonymous function in terms of both the dependent variable(s) and time:

>> func = @(t,y) y;
>> ode45(func, [0,1], 1)

Calling the ODE function without input arguments gives a graph of the solution. Calling it with one output argument returns a struct array:

>> Struct = ode45(func, [0,1], 1)
 Struct = solver: 'ode45'
          extdata: [1x1 struct]
          x: [1x11 double]
          y: [1x11 double]
          stats: [1x1 struct]
          idata: [1x1 struct]

This data is mostly used to, in the future, call the 'deval' function to get the answer at any time you want:

>> Solution = deval(Struct, 0.5)
 Solution = 1.6487

Deval can also return the derivative at the point of interest by including a second output argument.

>> [Solution, Derivative] = deval(Struct, 0.5)
Solution = 1.6487
Derivative = 1.6487

Since the derivative of e^x is itself it makes sense that the derivative and solution are the same here.

Calling ode45 with two output arguments returns two lists of data; t first, then the independent variables in an appropriately-sized matrix.

ODE Options

The are a rather large number of options that MATLAB gives you to modify how it solves the differential equations. The help file does a pretty good job describing all of them so they won't be described in detail here. To get a list use the help function:

>>help odeset

To get a list of the different options' names and what you have to pass to it, just type 'odeset' into the command prompt. It returns either a data type or a finite list of options. It also lists, in parenthesis, the default values of all the options.

To set a specific option or list of options, type the name of the option first and then the value of the option you want. For example, suppose you want to tighten the default relative tolerance from 10^-3 to 10^-4. You would call 'odeset' as follows:

>> option = odeset('RelTol', 10^-4);

Note that the option name must be passed as a string, or else you'll get an 'undefined variable' error most likely. More than one option can be passed at a time by just putting them all in a list:

>> option = odeset('RelTol', 10^-4, 'AbsTol', 10^-7, 'MassSingular', 'no');

The options structure can then be passed to the ode function as a forth (optional) input argument:

>> [T,X] = ode45(@func, [0,1], 1, option);

This will return more accurate values than default because the error tolerances are tighter. It should also compute faster because MATLAB is not checking to see if this is a differential-algebraic equation (this is what the MassSinglular option does; it is usually set to 'maybe' so MATLAB checks by itself).

The ODE solvers

MATLAB doesn't just have one ODE solver, it has eight as of the MATLAB 7.1 release. Some are more suited for certain problems than others, which is why all of them are included. The MATLAB help has a list of what functions each one can do, but here is a quick summary, in roughly the order you should try them unless you already know the classification of the problem:

Most problems can be solved using ode45, and since this is the best tradeoff of speed and accuracy it should be the first one you use.

If you need a really tight error tolerance or a lot of data points, use ode113.

If you have a relatively loose error tolerance or the problem is slow with ode45, try ode23.

If the problem is truly stiff and ode45 fails, use ode23tb for loose tolerances, ode23s for slightly tighter tolerances with a constant mass matrix, or ode15s for tighter tolerances or nonconstant mass matrices.

If you need a solution without numerical damping on a stiff problem, use ode23t.

The s indicates that the algorithm is intended for stiff problems.

There is one other ODE solver, which is special:

Implicit problems can only be solved using ode15i.

Since ode15i has some differences in its syntax, it is discussed next.

Implicit ODEs: ODE15i

Since ode15i is the only ODE solver that solves implicit equations, it must have some special syntactical rules on how to input the function.

If you are using ode15i declare the function as follows:

function Z = func(t,y,Yprime)
Z = 0;
Z = y - Yprime;

Note that with ode15i, you must put the function into normal form (solve it for 0), whereas for all other ODE functions you must solve explicitly for y'. Also notice that you must declare three input arguments instead of the usual two.

When you call ode15i, you must not only include initial conditions for y but also for Yprime. The initial conditions for Yprime go into the fourth argument:

>> [t,x] = ode15i(@func, [0,1], 1, 1);

This will return similar results to ode45 when used for the explicit equation, but has fewer data points.

The options structure is passed to ode15i as the optional fifth argument. All output options from ode15i are the same as for the other ODE solvers.