The following table lists the initial value problem solvers, the kind of problem you can solve with each solver, and the method each solver uses.
Solves These Kinds of Problems
Nonstiff differential equations
Nonstiff differential equations
Nonstiff differential equations
Stiff differential equations and DAEs
Stiff differential equations
Moderately stiff differential equations and DAEs
Stiff differential equations
Fully implicit differential equations
You can use the following functions to evaluate and extend solutions to ODEs.
Evaluate the numerical solution using the output of ODE solvers.
Extend the solution of an initial value problem for an ODE
An options structure contains named properties whose values are passed to ODE solvers, and which affect problem solution. Use these functions to create, alter, or access an options structure.
Create or alter options structure for input to ODE solver.
Extract properties from options structure created with
If an output function is specified, the solver calls the specified
function after every successful integration step. You can use
odeset to specify one of these sample
functions as the OutputFcn property, or
you can modify them to create your own functions.
An ordinary differential equation (ODE) contains one or more derivatives of a dependent variable y with respect to a single independent variable t, usually referred to as time. The derivative of y with respect to t is denoted as y ′, the second derivative as y ′′, and so on. Often y(t) is a vector, having elements y1, y2, ..., yn.
MATLAB® solvers handle the following types of first-order ODEs:
Explicit ODEs of the form y ′ = f (t, y)
Linearly implicit ODEs of the form M(t, y) y ′ = f (t, y), where M(t, y) is a matrix
Fully implicit ODEs of the form f (t, y, y′)
= 0 (
MATLAB ODE solvers accept only first-order differential equations. To use the solvers with higher-order ODEs, you must rewrite each equation as an equivalent system of first-order differential equations of the form
y′ = f(t,y)
You can write any ordinary differential equation
y(n) = f(t,y,y′,...,y(n − 1))
as a system of first-order equations by making the substitutions
y1 = y, y2 = y′,..., yn = y(n − 1)
y1= y, y2 = y', ... , yn = y(n − 1)
The result is an equivalent system of n first-order ODEs.
For example, you can rewrite the second-order van der Pol equation
by making the substitution
The resulting system of first-order ODEs is
Generally there are many functions y(t) that satisfy a given ODE, and additional information is necessary to specify the solution of interest. In an initial value problem, the solution of interest satisfies a specific initial condition, that is, y is equal to y0 at a given initial time t0. An initial value problem for an ODE is then
If the function f (t, y) is sufficiently smooth, this problem has one and only one solution. Generally there is no analytic expression for the solution, so it is necessary to approximate y(t) by numerical means.
There are three solvers designed for nonstiff problems:
Based on an
explicit Runge-Kutta (4,5) formula, the Dormand-Prince pair. It is
a one-step solver – in computing y(tn),
it needs only the solution at the immediately preceding time point, y(tn–1).
Based on an
explicit Runge-Kutta (2,3) pair of Bogacki and Shampine. It may be
more efficient than
Variable order Adams-Bashforth-Moulton PECE solver.
It may be more efficient than
Not all difficult problems are stiff, but all stiff problems are difficult for solvers not specifically designed for them. Solvers for stiff problems can be used exactly like the other solvers. However, you can often significantly improve the efficiency of these solvers by providing them with additional information about the problem. (See Integrator Options.)
There are four solvers designed for stiff problems:
solver based on the numerical differentiation formulas (NDFs). Optionally
it uses the backward differentiation formulas, BDFs (also known as
Gear's method). Like
on a modified Rosenbrock formula of order 2. Because it is a one-step
solver, it may be more efficient than
An implementation of the trapezoidal rule using a "free" interpolant. Use this solver if the problem is only moderately stiff and you need a solution without numerical damping.
An implementation of TR-BDF2, an implicit Runge-Kutta
formula with a first stage that is a trapezoidal rule step and a second
stage that is a backward differentiation formula of order 2. Like
fully implicit differential equations of the form
f(t,y,y') = 0
using the variable order BDF method. The basic syntax for
[t,y] = ode15i(odefun,tspan,y0,yp0,options)
The input arguments are
A function that evaluates the left side of the differential equation of the form f(t,y,y') = 0.
A vector specifying the interval of integration,
Vectors of initial conditions for y(t0)
respectively. The specified values must be consistent; that is, they
Optional integration argument created using the
The output arguments contain the solution approximated at discrete points:
Column vector of time points
Solution array. Each row in
page for more information about these arguments.
All of the ODE solver functions, except for
share a syntax that makes it easy to try any of the different numerical
methods, if it is not apparent which is the most appropriate. To apply
a different method to the same problem, simply change the ODE solver
function name. The simplest syntax, common to all the solver functions,
[t,y] = solver(odefun,tspan,y0,options)
solver is one of the ODE solver
functions listed previously.
The basic input arguments are
function handle that evaluates the system of ODEs. The function has the form
dydt = odefun(t,y)
Vector specifying the interval of integration. The solver
imposes the initial conditions at
Vector of initial conditions for the problem
See also Initial Value Problems.
Structure of optional parameters that change the default integration properties.
Integrator Options tells you how to create the structure and describes the properties you can specify.
The output arguments contain the solution approximated at discrete points:
Column vector of time points
Solution array. Each row in
See the reference page for the ODE solvers for more information about these arguments.
The default integration properties in the ODE solvers are selected
to handle common problems. In some cases, you can improve ODE solver
performance by overriding these defaults. You do this by supplying
the solvers with an
options structure that specifies
one or more property values.
For example, to change the value of the relative error tolerance
of the solver from the default value of
Create an options structure using the function
odeset by entering
options = odeset('RelTol', 1e-4);
to the solver as follows:
For all solvers except
use the syntax
[t,y] = solver(odefun,tspan,y0,options)
ode15i, use the syntax
[t,y] = ode15i(odefun,tspan,y0,yp0,options)
For a complete description of the available options, see the
reference page for
This example illustrates the steps for solving an initial value ODE problem:
Rewrite the problem as a system of first-order ODEs. Rewrite the van der Pol equation (second-order)
where μ > 0 is a scalar parameter, by making the substitution y'1 = y2. The resulting system of first-order ODEs is
Code the system of first-order ODEs. Once you represent the equation as a system of first-order ODEs, you can code it as a function that an ODE solver can use. The function must be of the form
dydt = odefun(t,y)
be the function's two arguments, the function does not need to use
them. The output
dydt, a column vector, is the
The code below represents the van der Pol system in the function,
vdp1 function assumes that μ = 1. The variables y1 and y2 are
a two-element vector.
function dydt = vdp1(t,y) dydt = [y(2); (1-y(1)^2)*y(2)-y(1)];
Note that, although
vdp1 must accept the
y, it does not
t in its computations.
Decide which solver you want to use to solve the problem. Then call the solver and pass it the function you created to describe the first-order system of ODEs, the time interval on which you want to solve the problem, and an initial condition vector.
For the van der Pol system, you can use
[0 20] with initial values
= 2 and
y(2) = 0.
[t,y] = ode45(@vdp1,[0 20],[2; 0]);
This example uses
vdp1 as a function handle to
The resulting output is a column vector of time points
a solution array
y. Each row in
to a time returned in the corresponding row of
The first column of
y corresponds to y1,
and the second column to y2.
View the solver output. You
can simply use the
to view the solver output.
plot(t,y(:,1),'-',t,y(:,2),'--') title('Solution of van der Pol Equation, \mu = 1'); xlabel('time t'); ylabel('solution y'); legend('y_1','y_2')
As an alternative, you can use a solver output function to process
the output. The solver calls the function specified in the integration
OutputFcn after each successful time step.
odeset to set
the desired function. See Solver Output
Properties, in the reference page for
for more information about
This example presents a stiff problem. For a stiff problem, solutions can change on a time scale that is very short compared to the interval of integration, but the solution of interest changes on a much longer time scale. Methods not designed for stiff problems are ineffective on intervals where the solution changes slowly because they use time steps small enough to resolve the fastest possible change.
When μ is increased to 1000, the solution
to the van der Pol equation changes dramatically and exhibits oscillation
on a much longer time scale. Approximating the solution of the initial
value problem becomes a more difficult task. Because this particular
problem is stiff, a solver intended for nonstiff problems, such as
ode45, is too inefficient to be practical.
A solver such as
intended for such stiff problems.
vdp1000 function evaluates the van
der Pol system from the previous example, but with μ = 1000.
function dydt = vdp1000(t,y) dydt = [y(2); 1000*(1-y(1)^2)*y(2)-y(1)];
Now use the
ode15s function to solve the
problem with the initial condition vector of
but a time interval of
[0 3000]. For scaling reasons,
plot just the first component of
[t,y] = ode15s(@vdp1000,[0 3000],[2; 0]); plot(t,y(:,1),'-'); title('Solution of van der Pol Equation, \mu = 1000'); xlabel('time t'); ylabel('solution y_1');
The preceding sections showed how to solve the van der Pol equation for two different values of the parameter µ. In those examples, the values µ = 1 and µ=1000 are hard-coded in the ODE functions. If you are solving an ODE for several different parameter values, it might be more convenient to include the parameter in the ODE function and assign a value to the parameter each time you run the ODE solver. This section explains how to do this for the van der Pol equation.
One way to provide parameter values to the ODE function is to write a MATLAB file that
Accepts the parameters as inputs.
Contains ODE function as a nested function, internally using the input parameters.
Calls the ODE solver.
The following code illustrates this:
function [t,y] = solve_vdp(mu) tspan = [0 max(20, 3*mu)]; y0 = [2; 0]; % Call the ODE solver ode15s. [t,y] = ode15s(@vdp,tspan,y0); % Define the ODE function as nested function, % using the parameter mu. function dydt = vdp(t,y) dydt = [y(2); mu*(1-y(1)^2)*y(2)-y(1)]; end end
Because the ODE function
vdp is a nested
function, the value of the parameter
mu is available
To run the MATLAB file for
mu = 1, enter
[t,y] = solve_vdp(1);
To run the code for µ = 1000, enter
[t,y] = solve_vdp(1000);
The numerical methods implemented in the ODE solvers produce
a continuous solution over the interval of integration [a,b].
You can evaluate the approximate solution, S(x),
at any point in [a,b] using
deval and the structure
by the solver. For example, if you solve the problem described in van der Pol Equation (Nonstiff) by calling
a single output argument
sol = ode45(@vdp1,[0 20],[2; 0]);
ode45 returns the solution as a structure.
You can then evaluate the approximate solution at points in the vector
= 1:5 as follows:
xint = 1:5; Sxint = deval(sol,xint) Sxint = 1.5081 0.3235 -1.8686 -1.7407 -0.8344 -0.7803 -1.8320 -1.0220 0.6260 1.3095
deval function is vectorized. For a
Sxint approximates the solution y(
rigidode illustrates the solution of a standard
test problem proposed by Krogh for solvers intended for nonstiff problems .
The ODEs are the Euler equations of a rigid body without external forces.
For your convenience, the entire problem is defined and solved
in a single MATLAB file. The differential equations are coded
as the local function
f. Because the example calls
ode45 solver without output arguments, the
solver uses the default output function
plot the solution components.
function rigidode %RIGIDODE Euler equations: rigid body without external forces tspan = [0 12]; y0 = [0; 1; 1]; % Solve the problem using ode45 ode45(@f,tspan,y0); % ------------------------------------------------------------ function dydt = f(t,y) dydt = [ y(2)*y(3) -y(1)*y(3) -0.51*y(1)*y(2) ];
The following example shows how to use the function
ode15i to solve the implicit ODE problem
defined by Weissinger's equation
ty2(y')3 – y3(y')2 + t(t2 + 1)y' – t2y = 0
with the initial value y(1) = (3/2)1/2. The exact solution of the ODE is
y(t) = (t2 + 0.5)1/2
The example uses the function
which is provided with MATLAB, to compute the left-hand side
of the equation.
ode15i, the example uses a
decic to compute a consistent initial
value for y'(t0).
In the following call, the given initial value y(1) = (3/2)1/2 is
held fixed and a guess of 0 is specified for y'(1).
See the reference page for
t0 = 1; y0 = sqrt(3/2); yp0 = 0; [y0,yp0] = decic(@weissinger,t0,y0,1,yp0,0);
You can now call
ode15i to solve the ODE
and then plot the numerical solution against the analytical solution
with the following commands.
[t,y] = ode15i(@weissinger,[1 10],y0,yp0); ytrue = sqrt(t.^2 + 0.5); plot(t,y,t,ytrue,'o');
fem1ode illustrates the solution of ODEs
that result from a finite element discretization of a partial differential
equation. The value of
N in the call
the discretization, and the resulting system consists of
N is 1
This example involves a mass matrix. The system of ODEs comes from a method of lines solution of the partial differential equation
with initial condition
u(0,x) = sin(x)
and boundary conditions
u(t,0) = u(t,π) = 0.
An integer N is chosen, h is defined as π/(N + 1), and the solution of the partial differential equation is approximated at xk = kh for k = 0, 1, ..., N+1 by
Here ϕk(x) is a piecewise linear function that is 1 atxk and 0 at all the otherxj. A Galerkin discretization leads to the system of ODEs
and the tridiagonal matrices M(t) and J are given by
The initial values c(0) are taken from the initial condition for the partial differential equation. The problem is solved on the time interval [0,π].
fem1ode example, the properties
options = odeset('Mass',@mass,'MStateDep','none','Jacobian',J)
indicate that the problem is of the form M(t)y' = Jy.
The nested function
mass(t) evaluates the time-dependent
mass matrix M(t) and
the constant Jacobian.
function fem1ode(N) %FEM1ODE Stiff problem with a time-dependent mass matrix if nargin < 1 N = 19; end h = pi/(N+1); y0 = sin(h*(1:N)'); tspan = [0; pi]; % The Jacobian is constant. e = repmat(1/h,N,1); % e=[(1/h) ... (1/h)]; d = repmat(-2/h,N,1); % d=[(-2/h) ... (-2/h)]; % J is shared with the derivative function. J = spdiags([e d e], -1:1, N, N); d = repmat(h/6,N,1); % M is shared with the mass matrix function. M = spdiags([d 4*d d], -1:1, N, N); options = odeset('Mass',@mass,'MStateDep','none', ... 'Jacobian',J); [t,y] = ode15s(@f,tspan,y0,options); figure; surf((1:N)/(N+1),t,y); set(gca,'ZLim',[0 1]); view(142.5,30); title(['Finite element problem with time-dependent mass ' ... 'matrix, solved by ODE15S']); xlabel('space ( x/\pi )'); ylabel('time'); zlabel('solution'); %-------------------------------------------------------------- function yp = f(t,y) % Derivative function. yp = J*y; % Constant Jacobian provided by outer function end % End nested function f %-------------------------------------------------------------- function Mt = mass(t) % Mass matrix function. Mt = exp(-t)*M; % M is provided by outer function end % End nested function mass %-------------------------------------------------------------- end
brussode illustrates the solution of a potentially
large stiff sparse problem. The problem is the classic "Brusselator"
system  that
models diffusion in a chemical reaction
and is solved on the time interval
[0,10] with α = 1/50 and
ui(0) = 1
vi(0) = 3
where, for i = 1,...,N, xi = i/(N + 1). There are 2N equations in the system, but the Jacobian is banded with a constant width 5 if the equations are ordered as u1, v1, u2, v2, ...
In the call
to N, the parameter
the number of grid points. The resulting system consists of
20. The problem
becomes increasingly stiff and the Jacobian increasingly sparse as
The nested function
f(t,y) returns the derivatives
vector for the Brusselator problem. The local function
a sparse matrix of
showing the locations of nonzeros in the Jacobian ∂f/∂y.
The example assigns this matrix to the property
and the solver uses the sparsity pattern to generate the Jacobian
numerically as a sparse matrix. Providing a sparsity pattern can significantly
reduce the number of function evaluations required to generate the
Jacobian and can accelerate integration.
For the Brusselator problem, if the sparsity pattern is not
N evaluations of the function are needed
to compute the 2
matrix. If the sparsity pattern is supplied, only four evaluations
are needed, regardless of the value of
function brussode(N) %BRUSSODE Stiff problem modeling a chemical reaction if nargin < 1 N = 20; end tspan = [0; 10]; y0 = [1+sin((2*pi/(N+1))*(1:N)); repmat(3,1,N)]; options = odeset('Vectorized','on','JPattern',jpattern(N)); [t,y] = ode15s(@f,tspan,y0,options); u = y(:,1:2:end); x = (1:N)/(N+1); surf(x,t,u); view(-40,30); xlabel('space'); ylabel('time'); zlabel('solution u'); title(['The Brusselator for N = ' num2str(N)]); % -------------------------------------------------------------- function dydt = f(t,y) c = 0.02 * (N+1)^2; dydt = zeros(2*N,size(y,2)); % preallocate dy/dt % Evaluate the two components of the function at one edge of % the grid (with edge conditions). i = 1; dydt(i,:) = 1 + y(i+1,:).*y(i,:).^2 - 4*y(i,:) + ... c*(1-2*y(i,:)+y(i+2,:)); dydt(i+1,:) = 3*y(i,:) - y(i+1,:).*y(i,:).^2 + ... c*(3-2*y(i+1,:)+y(i+3,:)); % Evaluate the two components of the function at all interior % grid points. i = 3:2:2*N-3; dydt(i,:) = 1 + y(i+1,:).*y(i,:).^2 - 4*y(i,:) + ... c*(y(i-2,:)-2*y(i,:)+y(i+2,:)); dydt(i+1,:) = 3*y(i,:) - y(i+1,:).*y(i,:).^2 + ... c*(y(i-1,:)-2*y(i+1,:)+y(i+3,:)); % Evaluate the two components of the function at the other edge % of the grid (with edge conditions). i = 2*N-1; dydt(i,:) = 1 + y(i+1,:).*y(i,:).^2 - 4*y(i,:) + ... c*(y(i-2,:)-2*y(i,:)+1); dydt(i+1,:) = 3*y(i,:) - y(i+1,:).*y(i,:).^2 + ... c*(y(i-1,:)-2*y(i+1,:)+3); end % End nested function f end % End function brussode % -------------------------------------------------------------- function S = jpattern(N) B = ones(2*N,5); B(2:2:2*N,2) = zeros(N,1); B(1:2:2*N-1,4) = zeros(N,1); S = spdiags(B,-2:2,2*N,2*N); end;
ballode models the motion of a bouncing ball.
This example illustrates the event location capabilities of the ODE
The equations for the bouncing ball are:
y'1 = y2
y'2 = –9.8
In this example, the event function is coded in the local function
[value,isterminal,direction] = events(t,y)
A value of the event function
The information whether or not the integration should
value = 0 (
isterminal = 1 or
The desired directionality of the zero crossings:
Detect zero crossings in the negative direction only
Detect all zero crossings
Detect zero crossings in the positive direction only
The length of
direction is the same as the number of event
ith element of each vector, corresponds
ith event function. For an example of more
advanced event location, see
orbitode (Advanced Event Location).
ballode, setting the
@events causes the solver to stop the integration
isterminal = 1) when the ball hits the ground
its fall (
direction = -1). The example then restarts
the integration with initial conditions corresponding to a ball that
function ballode tstart = 0; tfinal = 30; y0 = [0; 20]; refine = 4; options = odeset('Events',@events,'OutputFcn', @odeplot,... 'OutputSel',1,'Refine',refine); set(gca,'xlim',[0 30],'ylim',[0 25]); box on hold on; tout = tstart; yout = y0.'; teout = ; yeout = ; ieout = ; for i = 1:10 % Solve until the first terminal event. [t,y,te,ye,ie] = ode23(@f,[tstart tfinal],y0,options); if ~ishold hold on end % Accumulate output. nt = length(t); tout = [tout; t(2:nt)]; yout = [yout; y(2:nt,:)]; teout = [teout; te]; % Events at tstart are never reported. yeout = [yeout; ye]; ieout = [ieout; ie]; ud = get(gcf,'UserData'); if ud.stop break; end % Set the new initial conditions, with .9 attenuation. y0(1) = 0; y0(2) = -.9*y(nt,2); % A good guess of a valid first time step is the length of % the last valid time step, so use it for faster computation. options = odeset(options,'InitialStep',t(nt)-t(nt-refine),... 'MaxStep',t(nt)-t(1)); tstart = t(nt); end plot(teout,yeout(:,1),'ro') xlabel('time'); ylabel('height'); title('Ball trajectory and the events'); hold off odeplot(,,'done'); % -------------------------------------------------------------- function dydt = f(t,y) dydt = [y(2); -9.8]; % -------------------------------------------------------------- function [value,isterminal,direction] = events(t,y) % Locate the time when height passes through zero in a % decreasing direction and stop integration. value = y(1); % Detect height = 0 isterminal = 1; % Stop the integration direction = -1; % Negative direction only
orbitode illustrates the solution of a standard
test problem for those solvers that are intended for nonstiff problems.
It traces the path of a spaceship traveling around the moon and returning
to the earth (Shampine and Gordon , p. 246).
orbitode problem is a system of the following
The first two solution components are coordinates of the body
of infinitesimal mass, so plotting one against the other gives the
orbit of the body. The initial conditions have been chosen to make
the orbit periodic. The value of μ corresponds
to a spaceship traveling around the moon and the earth. Moderately
stringent tolerances are necessary to reproduce the qualitative behavior
of the orbit. Suitable values are
events function includes event
functions that locate the point of maximum distance from the starting
point and the time the spaceship returns to the starting point. Note
that the events are located accurately, even though the step sizes
used by the integrator are not determined by
the location of the events. In this example, the ability to specify
the direction of the zero crossing is critical. Both the point of
return to the initial point and the point of maximum distance have
the same event function value, and the direction of the crossing is
used to distinguish them.
To run this example, type
the command line. The example uses the output function
produce the two-dimensional phase plane plot and let you to see the
progress of the integration.
function orbitode %ORBITODE Restricted three-body problem mu = 1 / 82.45; mustar = 1 - mu; y0 = [1.2; 0; 0; -1.04935750983031990726]; tspan = [0 7]; options = odeset('RelTol',1e-5,'AbsTol',1e-4,... 'OutputFcn',@odephas2,'Events',@events); [t,y,te,ye,ie] = ode45(@f,tspan,y0,options); plot(y(:,1),y(:,2),ye(:,1),ye(:,2),'o'); title ('Restricted three body problem') ylabel ('y(t)') xlabel ('x(t)') % -------------------------------------------------------------- function dydt = f(t,y) r13 = ((y(1) + mu)^2 + y(2)^2) ^ 1.5; r23 = ((y(1) - mustar)^2 + y(2)^2) ^ 1.5; dydt = [ y(3) y(4) 2*y(4) + y(1) - mustar*((y(1)+mu)/r13) - ... mu*((y(1)-mustar)/r23) -2*y(3) + y(2) - mustar*(y(2)/r13) - mu*(y(2)/r23) ]; end % End nested function f % -------------------------------------------------------------- function [value,isterminal,direction] = events(t,y) % Locate the time when the object returns closest to the % initial point y0 and starts to move away; stop integration. % Also locate the time when the object is farthest from the % initial point y0 and starts to move closer. % % The current distance of the body is % % DSQ = (y(1)-y0(1))^2 + (y(2)-y0(2))^2 % = <y(1:2)-y0(1:2),y(1:2)-y0(1:2)> % % A local minimum of DSQ occurs when d/dt DSQ crosses zero % heading in the positive direction. Compute d(DSQ)/dt as % % d(DSQ)/dt = 2*(y(1:2)-y0(1:2))'*dy(1:2)/dt = ... % 2*(y(1:2)-y0(1:2))'*y(3:4) % dDSQdt = 2 * ((y(1:2)-y0(1:2))' * y(3:4)); value = [dDSQdt; dDSQdt]; isterminal = [1; 0]; % Stop at local minimum direction = [1; -1]; % [local minimum, local maximum] end % End nested function events end
Note: The Robertson problem appears as an example in the prolog to LSODI .
hb1ode, the problem is solved with initial
conditions y1(0) = 1, y2(0) = 0, y3(0) = 0 to steady state.
These differential equations satisfy a linear conservation law that
is used to reformulate the problem as the DAE
These equations do not have a solution for y(0) with components that do not sum to 1. The problem has the form of My' = f(t,y) with
M is singular, but
not inform the solver of this. The solver must recognize that the
problem is a DAE, not an ODE. Similarly, although consistent initial
conditions are obvious, the example uses an inconsistent value y3(0) = 10–3 to
illustrate computation of consistent initial conditions.
Imposes a much smaller absolute error tolerance on y2 than on the other components. This is because y2 is much smaller than the other components and its major change takes place in a relatively short time.
Specifies additional points at which the solution is computed to more clearly show the behavior of y2.
Multiplies y2 by 104 to make y2 visible when plotting it with the rest of the solution.
Uses a logarithmic scale to plot the solution on the long time interval.
function hb1dae %HB1DAE Stiff differential-algebraic equation (DAE) % A constant, singular mass matrix M = [1 0 0 0 1 0 0 0 0]; % Use inconsistent initial condition to test initialization. y0 = [1; 0; 1e-3]; tspan = [0 4*logspace(-6,6)]; % Use the LSODI example tolerances.'MassSingular' is left % at its default 'maybe' to test the automatic detection % of a DAE. options = odeset('Mass',M,'RelTol',1e-4,... 'AbsTol',[1e-6 1e-10 1e-6],... 'Vectorized','on'); [t,y] = ode15s(@f,tspan,y0,options); y(:,2) = 1e4*y(:,2); semilogx(t,y); ylabel('1e4 * y(:,2)'); title(['Robertson DAE problem with a Conservation Law, '... 'solved by ODE15S']); xlabel('This is equivalent to the stiff ODEs in HB1ODE.'); % -------------------------------------------------------- function out = f(t,y) out = [ -0.04*y(1,:) + 1e4*y(2,:).*y(3,:) 0.04*y(1,:) - 1e4*y(2,:).*y(3,:) - 3e7*y(2,:).^2 y(1,:) + y(2,:) + y(3,:) - 1 ];
If certain components of the solution must be nonnegative, use
odeset to set the
for the indices of these components.
Imposing nonnegativity is not always a trivial task. We suggest that you use this option only when necessary, for example in instances in which the application of a solution or integration will fail otherwise.
Consider the following initial value problem solved on the interval
y' = - |y|, y(0) = 1
The solution of this problem decays to zero. If a solver produces
a negative approximate solution, it begins to track the solution of
the ODE through this value, the solution goes off to minus infinity,
and the computation fails. Using the
prevents this from happening.
In this example, the first call to
the defaults for the solver parameters:
ode = @(t,y) -abs(y); [t0,y0] = ode45(ode,[0, 40], 1);
The second uses options to impose nonnegativity conditions:
options = odeset('NonNegative',1); [t1,y1] = ode45(ode,[0, 40], 1, options);
This plot compares the numerical solution to the exact solution.
Here is a more complete view of the code used to obtain this plot:
ode = @(t,y) -abs(y); options = odeset('Refine',1); [t0,y0] = ode45(ode,[0, 40], 1,options); options = odeset(options,'NonNegative',1); [t1,y1] = ode45(ode,[0, 40], 1, options); t = linspace(0,40,1000); y = exp(-t); plot(t,y,'b-',t0,y0,'ro',t1,y1,'b*'); legend('Exact solution','No constraints','Nonnegativity', ... 'Location','SouthWest')
The kneeode Example. The
kneeode example solves the "knee
problem" by imposing a nonnegativity constraint on the numerical
solution. The initial value problem is
ε*y' = (1-x)*y - y^2, y(0) = 1
0 < ε
the solution of this problem approaches null isoclines
= 1 - x and
y = 0 for
x > 1, respectively. The numerical
solution, when computed with default tolerances, follows the
= 1 - x isocline for the whole interval of integration.
Imposing nonnegativity constraints results in the correct solution.
Here is the code that makes up the
function kneeode %KNEEODE The "knee problem" with Nonnegativity constraints. % Problem parameter epsilon = 1e-6; y0 = 1; xspan = [0, 2]; % Solve without imposing constraints options = ; [x1,y1] = ode15s(@odefcn,xspan,y0,options); % Impose nonnegativity constraint options = odeset('NonNegative',1); [x2,y2] = ode15s(@odefcn,xspan,y0,options); figure plot(x1,y1,'b.-',x2,y2,'g-') axis([0,2,-1,1]); title('The "knee problem"'); legend('No constraints','nonnegativity') xlabel('x'); ylabel('solution y') function yp = odefcn(x,y) yp = ((1 - x)*y - y^2)/epsilon; end end % kneeode
The derivative function is defined within nested function
The value of epsilon used in
odefcn is obtained
from the outer function:
function yp = odefcn(x,y) yp = ((1 - x)*y - y^2)/epsilon; end
The example solves the problem using the
ode15s function, first with the default
options, and then by imposing a nonnegativity constraint. To run the
kneeode at the MATLAB command
Here is the output plot. The plot confirms correct solution behavior after imposing constraints.
The following additional examples are available. Type
to view the code and
to run the example.
Stiff DAE — electrical circuit
Simple event location — bouncing ball
ODE with time- and state-dependent mass matrix — motion of a baton
Stiff large problem — diffusion in a chemical reaction (the Brusselator)
ODE with strongly state-dependent mass matrix — Burgers' equation solved using a moving mesh technique
Stiff problem with a time-dependent mass matrix — finite element method
Stiff problem with a constant mass matrix — finite element method
Stiff ODE problem solved on a very long interval — Robertson chemical reaction
Robertson problem — stiff, linearly implicit DAE from a conservation law
Robertson problem — stiff, fully implicit DAE
Burgers' equation solved as implicit ODE system
The "knee problem" with nonnegativity constraints
Advanced event location — restricted three body problem
Nonstiff problem — Euler equations of a rigid body without external forces
Parameterizable van der Pol equation (stiff for large μ)
do the ODE solvers differ from
Can I solve ODE systems in which there are more equations than unknowns, or vice versa?
Memory and Computational Efficiency
How large a problem can I solve with the ODE suite?
The primary constraints are memory and time. At each
time step, the solvers for nonstiff problems allocate vectors of length
If the problem is nonstiff, or if you are using the sparse option, it may be possible to solve a problem with thousands of unknowns. In this case, however, storage of the result can be problematic. Try asking the solver to evaluate the solution at specific points only, or call the solver with no output arguments and use an output function to monitor the solution.
solving a very large system, but only care about a couple of the components
Yes. The user-installable
output function capability is designed specifically for this
purpose. When you call the solver with no output arguments, the solver
does not allocate storage to hold the entire solution history. Instead,
the solver calls
What is the startup cost of the integration and how can I reduce it?
The biggest startup cost occurs as the solver attempts
to find a step size appropriate to the scale of the problem. If you
happen to know an appropriate step size, use the
Time Steps for Integration
The first step size that the integrator takes is too large, and it misses important behavior.
You can specify the first step size with the
Can I integrate with fixed step sizes?
Error Tolerance and Other Options
speaking, this means that you want
I want answers that are correct
to the precision of the computer. Why can't I simply set
You can get close to machine precision, but not that
close. The solvers do not allow
How do I tell the solver that I don't care about getting an accurate answer for one of the solution components?
You can increase the absolute error tolerance corresponding to this solution component. If the tolerance is bigger than the component, this specifies no correct digits for the component. The solver may have to get some correct digits in this component to compute other components accurately, but it generally handles this automatically.
Other Types of Equations
Can the solvers handle partial differential equations (PDEs) that have been discretized by the method of lines?
Yes, because the discretization produces a system
of ODEs. Depending on the discretization, you might have a form involving
mass matrices – the ODE solvers provide for this. Often the
system is stiff. This is to be expected when the PDE is parabolic
and when there are phenomena that happen on very different time scales
such as a chemical reaction in a fluid flow. In such cases, use one
of the four solvers:
there are many equations, set the
Can I solve differential-algebraic equation (DAE) systems?
Yes. The solvers
Can I integrate a set of sampled data?
Not directly. You have to represent the data as a function
by interpolation or some other scheme for fitting data. The smoothness
of this function is critical. A piecewise polynomial fit like a spline
can look smooth to the eye, but rough to a solver; the solver takes
small steps where the derivatives of the fit have jumps. Either use
a smooth function to represent the data or use one of the lower order
What do I do when I have the final and not the initial value?
All the solvers of the ODE suite allow you to solve backwards
or forwards in time. The syntax for the solvers is
Other Common Problems
The solution doesn't look like what I expected.
If you're right about its appearance, you need to reduce the error tolerances from their default values. A smaller relative error tolerance is needed to compute accurately the solution of problems integrated over "long" intervals, as well as solutions of problems that are moderately unstable.
You should check whether there are solution components that stay smaller than their absolute error tolerance for some time. If so, you are not asking for any correct digits in these components. This may be acceptable for these components, but failing to compute them accurately may degrade the accuracy of other components that depend on them.
My plots aren't smooth enough.
Increase the value of
I'm plotting the solution as it is computed and it looks fine, but the code gets stuck at some point.
First verify that the ODE function is smooth near the
point where the code gets stuck. If it isn't, the solver must take
small steps to deal with this. It may help to break
If the function
is smooth and the code is taking extremely small steps, you are probably
trying to solve a stiff problem with a solver not intended for this
purpose. Switch to
My integration proceeds very slowly, using too many time steps.
First, check that your
Finally, make sure that the ODE function is written in an efficient way. The solvers evaluate the derivatives in the ODE function many times. The cost of numerical integration depends critically on the expense of evaluating the ODE function. Rather than recompute complicated constant parameters at each evaluation, store them in globals or calculate them once and pass them to nested functions.
I know that the solution undergoes a radical change at time t where
t0 ≤ t ≤ tf
but the integrator steps past without "seeing" it.
If you know there is a sharp change at time
If the differential equation has periodic coefficients or solution, you might restrict the maximum step size to the length of the period so the integrator won't step over periods.