## Documentation Center |

Solve stiff differential equations and DAEs; variable order method

`[T,Y] = solver(odefun,tspan,y0)[T,Y] = solver(odefun,tspan,y0,options)[T,Y,TE,YE,IE] = solver(odefun,tspan,y0,options)sol = solver(odefun,[t0
tf],y0...)`

This page contains an overview of the solver functions: `ode23`, `ode45`, `ode113`, `ode15s`, `ode23s`, `ode23t`,
and `ode23tb`. You can call any of these solvers
by substituting the placeholder, * solver*,
with any of the function names.

The following table describes the input arguments to the solvers.

A function handle that evaluates the right side of the
differential equations. All solvers solve systems of equations in
the form | |

A vector specifying the interval of integration, For | |

Specifying Specifying | |

A vector of initial conditions. | |

Structure of optional parameters that change the default integration properties. This is the fourth input argument. [t,y] = You can create options using the |

The following table lists the output arguments for the solvers.

| Column vector of time points. |

| Solution array. Each row in |

| The time at which an event occurs. |

| The solution at the time of the event. |

| The index |

| Structure to evaluate the solution. |

`[T,Y] = solver(odefun,tspan,y0)`
with

Parameterizing Functions explains how to provide additional
parameters to the function `fun`, if necessary.

`[T,Y] = solver(odefun,tspan,y0,options)` solves
as above with default integration parameters replaced by property values specified
in

`[T,Y,TE,YE,IE] = solver(odefun,tspan,y0,options)` solves
as above while also finding where functions of (

`value(i)`is the value of the function.`isterminal(i) = 1`, if the integration is to terminate at a zero of this event function and`0`otherwise.`direction(i) = 0`if all zeros are to be computed (the default),`+1`if only the zeros where the event function increases, and`-1`if only the zeros where the event function decreases.

Corresponding entries in `TE`, `YE`,
and `IE` return, respectively, the time at which
an event occurs, the solution at the time of the event, and the index `i` of
the event function that vanishes.

`sol = solver(odefun,[t0
tf],y0...)` returns a structure that you can use with

Steps chosen by the solver. | |

Each column | |

Solver name. |

If you specify the `Events` option and events
are detected, `sol` also includes these fields:

Points at which events, if any, occurred. | |

Solutions that correspond to events in | |

Indices into the vector returned by the function specified
in the |

If you specify an output function as the value of the `OutputFcn` property,
the solver calls it with the computed solution after each time step.
Four output functions are provided: `odeplot`, `odephas2`, `odephas3`, `odeprint`.
When you call the solver with no output arguments, it calls the default `odeplot` to
plot the solution as it is computed. `odephas2` and `odephas3` produce
two- and three-dimensional phase plane plots, respectively. `odeprint` displays
the solution components on the screen. By default, the ODE solver
passes all components of the solution to the output function. You
can pass only specific components by providing a vector of indices
as the value of the `OutputSel` property. For example,
if you call the solver with no output arguments and set the value
of `OutputSel` to `[1,3]`, the solver
plots solution components 1 and 3 as they are computed.

For the stiff solvers `ode15s`, `ode23s`, `ode23t`,
and `ode23tb`, the Jacobian matrix ∂*f*/∂*y* is
critical to reliability and efficiency. Use `odeset` to
set `Jacobian` to `@FJAC` if `FJAC(T,Y)` returns
the Jacobian ∂*f*/∂*y* or
to the matrix ∂*f*/∂*y* if
the Jacobian is constant. If the `Jacobian` property
is not set (the default), ∂*f*/∂*y* is
approximated by finite differences. Set the `Vectorized` property
'`on`' if the ODE function is coded so that `odefun`(`T`,[`Y1,Y2
...`]) returns [`odefun`(`T`,`Y1`)`,odefun`(`T`,`Y2`) `...`].
If ∂*f*/∂*y* is a
sparse matrix, set the `JPattern` property to the
sparsity pattern of ∂*f*/∂*y*,
i.e., a sparse matrix `S` with `S(i,j) =
1` if the `i`th component of *f*(*t*,*y*)
depends on the `j`th component of *y*,
and 0 otherwise.

The solvers of the ODE suite can solve problems of the form *M*(*t*,*y*)*y*′ = *f*(*t*,*y*),
with time- and state-dependent mass matrix *M*. (The `ode23s` solver
can solve only equations with constant mass matrices.) If a problem
has a mass matrix, create a function `M = MASS(t,y)` that
returns the value of the mass matrix, and use `odeset` to
set the `Mass` property to `@MASS`.
If the mass matrix is constant, the matrix should be used as the value
of the `Mass` property. Problems with state-dependent
mass matrices are more difficult:

If the mass matrix does not depend on the state variable

*y*and the function`MASS`is to be called with one input argument,`t`, set the`MStateDependence`property to '`none`'.If the mass matrix depends weakly on

*y*, set`MStateDependence`to '`weak`' (the default); otherwise, set it to '`strong`'. In either case, the function`MASS`is called with the two arguments (`t`,`y`).

If there are many differential equations, it is important to exploit sparsity:

Return a sparse

*M*(*t*,*y*).Supply the sparsity pattern of ∂

*f*/∂*y*using the`JPattern`property or a sparse ∂*f*/∂*y*using the`Jacobian`property.For strongly state-dependent

*M*(*t*,*y*), set`MvPattern`to a sparse matrix`S`with`S(i,j) = 1`if for any`k`, the (`i,k`) component of*M*(*t*,*y*) depends on component`j`of*y*, and`0`otherwise.

If the mass matrix *M* is singular, then *M*(*t*,*y*)*y*′ = *f*(*t*,*y*)
is a system of differential algebraic equations. DAEs have solutions
only when *y*_{0} is consistent,
that is, if there is a vector *yp*_{0} such
that *M*(*t*_{0},*y*_{0})*yp*_{0} = *f*(*t*_{0},*y*_{0}).
The `ode15s` and `ode23t` solvers
can solve DAEs of index 1 provided that `y0` is sufficiently
close to being consistent. If there is a mass matrix, you can use `odeset` to set the `MassSingular` property
to `'yes'`, `'no'`, or `'maybe'`.
The default value of `'maybe'` causes the solver
to test whether the problem is a DAE. You can provide `yp0` as
the value of the `InitialSlope` property. The default
is the zero vector. If a problem is a DAE, and `y0` and `yp0` are
not consistent, the solver treats them as guesses, attempts to compute
consistent values that are close to the guesses, and continues to
solve the problem. When solving DAEs, it is very advantageous to formulate
the problem so that *M* is a diagonal matrix (a semi-explicit
DAE).

Solver | Problem Type | Order of Accuracy | When to Use |
---|---|---|---|

Nonstiff | Medium | Most of the time. This should be the first solver you try. | |

Nonstiff | Low | For problems with crude error tolerances or for solving moderately stiff problems. | |

Nonstiff | Low to high | For problems with stringent error tolerances or for solving computationally intensive problems. | |

Stiff | Low to medium | If | |

Stiff | Low | If using crude error tolerances to solve stiff systems and the mass matrix is constant. | |

Moderately Stiff | Low | For moderately stiff problems if you need a solution without numerical damping. | |

Stiff | Low | If using crude error tolerances to solve stiff systems. |

The algorithms used in the ODE solvers vary according to order of accuracy [6] and the type of systems (stiff or nonstiff) they are designed to solve. See Algorithms for more details.

Different solvers accept different parameters in the options
list. For more information, see `odeset` and Integrator Options in the MATLAB^{®} Mathematics
documentation.

Parameters | ode45 | ode23 | ode113 | ode15s | ode23s | ode23t | ode23tb |
---|---|---|---|---|---|---|---|

√ | √ | √ | √ | √ | √ | √ | |

√ | √ | √ | √ | √ | √ | √ | |

√ | √ | √ | √ * | — | √ * | √ * | |

√ | √ | √ | √ | √ | √ | √ | |

√ | √ | √ | √ | √ | √ | √ | |

— | — | — | √ | √ | √ | √ | |

√ | √ | √ | √ | √ | √ | √ | |

— | — | — | √ | — | √ | — | |

— | — | — | √ | — | — | — |

An example of a nonstiff system is the system of equations describing the motion of a rigid body without external forces.

To simulate this system, create a function `rigid` containing
the equations

function dy = rigid(t,y) dy = zeros(3,1); % a column vector dy(1) = y(2) * y(3); dy(2) = -y(1) * y(3); dy(3) = -0.51 * y(1) * y(2);

In this example we change the error tolerances using the `odeset` command and solve on a time interval `[0
12]` with an initial condition vector `[0 1 1]` at
time `0`.

options = odeset('RelTol',1e-4,'AbsTol',[1e-4 1e-4 1e-5]); [T,Y] = ode45(@rigid,[0 12],[0 1 1],options);

Plotting the columns of the returned array `Y` versus `T` shows
the solution

plot(T,Y(:,1),'-',T,Y(:,2),'-.',T,Y(:,3),'.')

An example of a stiff system is provided by the van der Pol equations in relaxation oscillation. The limit cycle has portions where the solution components change slowly and the problem is quite stiff, alternating with regions of very sharp change where it is not stiff.

To simulate this system, create a function `vdp1000` containing
the equations

function dy = vdp1000(t,y) dy = zeros(2,1); % a column vector dy(1) = y(2); dy(2) = 1000*(1 - y(1)^2)*y(2) - y(1);

For this problem, we will use the default relative and absolute
tolerances (`1e-3` and `1e-6`, respectively)
and solve on a time interval of `[0 3000]` with initial
condition vector `[2 0]` at time `0`.

[T,Y] = ode15s(@vdp1000,[0 3000],[2 0]);

Plotting the first column of the returned matrix `Y` versus `T` shows
the solution

plot(T,Y(:,1),'-o')

This example solves an ordinary differential equation with time-dependent terms.

Consider the following ODE, with time-dependent parameters defined only through the set of data points given in two vectors:

y'(t) + f(t)y(t) = g(t)

The
initial condition is `y(0) = 0`, where the function `f(t)` is
defined through the `n`-by-`1` vectors `tf` and `f`,
and the function `g(t)` is defined through the `m`-by-`1` vectors `tg` and `g`.

First, define the time-dependent parameters `f(t)` and `g(t)` as
the following:

ft = linspace(0,5,25); % Generate t for f f = ft.^2 - ft - 3; % Generate f(t) gt = linspace(1,6,25); % Generate t for g g = 3*sin(gt-0.25); % Generate g(t)

Write a function to interpolate the data sets specified above to obtain the value of the time-dependent terms at the specified time:

function dydt = myode(t,y,ft,f,gt,g) f = interp1(ft,f,t); % Interpolate the data set (ft,f) at time t g = interp1(gt,g,t); % Interpolate the data set (gt,g) at time t dydt = -f.*y + g; % Evaluate ODE at time t

Call the derivative function `myode.m` within
the MATLAB `ode45` function specifying time
as the first input argument :

Tspan = [1 5]; % Solve from t=1 to t=5 IC = 1; % y(t=0) = 1 [T Y] = ode45(@(t,y) myode(t,y,ft,f,gt,g),Tspan,IC); % Solve ODE

Plot the solution `y(t)` as a function
of time:

plot(T, Y); title('Plot of y as a function of time'); xlabel('Time'); ylabel('Y(t)');

[1] Bank, R. E., W. C. Coughran, Jr., W.
Fichtner, E. Grosse, D. Rose, and R. Smith, "Transient
Simulation of Silicon Devices and Circuits," *IEEE
Trans. CAD*, 4 (1985), pp. 436–451.

[2] Bogacki, P. and L. F. Shampine, "A
3(2) pair of Runge-Kutta formulas," *Appl. Math.
Letters*, Vol. 2, 1989, pp. 321–325.

[3] Dormand, J. R. and P. J. Prince, "A
family of embedded Runge-Kutta formulae," *J. Comp.
Appl. Math.*, Vol. 6, 1980, pp. 19–26.

[4] Forsythe, G. , M. Malcolm, and C. Moler, *Computer
Methods for Mathematical Computations*, Prentice-Hall,
New Jersey, 1977.

[5] Kahaner, D. , C. Moler, and S. Nash, *Numerical
Methods and Software*, Prentice-Hall, New Jersey, 1989.

[6] Shampine, L. F. , *Numerical
Solution of Ordinary Differential Equations*, Chapman
& Hall, New York, 1994.

[7] Shampine, L. F. and M. K. Gordon, *Computer
Solution of Ordinary Differential Equations: the Initial Value Problem*,
W. H. Freeman, San Francisco, 1975.

[8] Shampine, L. F. and M. E. Hosea, "Analysis
and Implementation of TR-BDF2," *Applied Numerical
Mathematics 20*, 1996.

[9] Shampine, L. F. and M. W. Reichelt, "The
MATLAB ODE Suite," *SIAM Journal on Scientific Computing, *Vol.
18, 1997, pp. 1–22.

[10] Shampine, L. F., M. W. Reichelt, and
J.A. Kierzenka, "Solving Index-1 DAEs in MATLAB and Simulink," *SIAM
Review*, Vol. 41, 1999, pp. 538–552.

`deval` | `function_handle` | `ode15i` | `odeget` | `odeset`

Was this topic helpful?