Solve stiff differential equations and DAEs — variable order method
tspan = [t0 tf], integrates the system of
differential equations from
y0. Each row in the solution
y corresponds to a value returned in column
All MATLAB® ODE solvers can solve systems of equations of
the form ,
or problems that involve a mass matrix, .
ode23s solver only can solve problems with
a mass matrix if the mass matrix is constant.
solve problems with a mass matrix that is singular, known as differential-algebraic
equations (DAEs). Specify the mass matrix using the
finds where functions of (t,y),
called event functions, are zero. In the output,
the time of the event,
ye is the solution at the
time of the event, and
ie is the index of the triggered
For each event function, specify whether the integration is
to terminate at a zero and whether the direction of the zero crossing
matters. Do this by setting the
to a function, such as
and creating a corresponding function: [
For more information, see ODE Event Location.
Simple ODEs that have a single solution component can be specified as an anonymous function in the call to the solver. The anonymous function must accept two inputs
(t,y) even if one of the inputs is not used.
Solve the ODE
Use a time interval of
[0,2] and the initial condition
y0 = 1.
tspan = [0 2]; y0 = 1; [t,y] = ode15s(@(t,y) -10*t, tspan, y0);
Plot the solution.
An example of a stiff system of equations is the van der Pol equations in relaxation oscillation. The limit cycle has regions where the solution components change slowly and the problem is quite stiff, alternating with regions of very sharp change where it is not stiff.
The system of equations is:
The initial conditions are
. The function
vdp1000 ships with MATLAB® and encodes the equations.
function dydt = vdp1000(t,y) %VDP1000 Evaluate the van der Pol ODEs for mu = 1000. % % See also ODE15S, ODE23S, ODE23T, ODE23TB. % Jacek Kierzenka and Lawrence F. Shampine % Copyright 1984-2014 The MathWorks, Inc. dydt = [y(2); 1000*(1-y(1)^2)*y(2)-y(1)];
Solving this system using
ode45 with the default relative and absolute error tolerances (
1e-6, respectively) is extremely slow, requiring several minutes to solve and plot the solution.
ode45 requires millions of time steps to complete the integration, due to the areas of stiffness where it struggles to meet the tolerances.
This is a plot of the solution obtained by
ode45, which takes a long time to compute. Notice the enormous number of time steps required to pass through areas of stiffness.
Solve the stiff system using the
ode15s solver, and then plot the first column of the solution
y against the time points
ode15s solver passes through stiff areas with far fewer steps than
[t,y] = ode15s(@vdp1000,[0 3000],[2 0]); plot(t,y(:,1),'-o')
ode15s only works with functions that use two input arguments,
y. However, you can pass in extra parameters by defining them outside the function and passing them in when you specify the function handle.
Solve the ODE
Rewriting the equation as a first-order system yields
odefcn.m represents this system of equations as a function that accepts four input arguments:
function dydt = odefcn(t,y,A,B) dydt = zeros(2,1); dydt(1) = y(2); dydt(2) = (A/B)*t.*y(1);
Solve the ODE using
ode15s. Specify the function handle such that it passes in the predefined values for
A = 1; B = 2; tspan = [0 5]; y0 = [0 0.01]; [t,y] = ode15s(@(t,y) odefcn(t,y,A,B), tspan, y0);
Plot the results.
ode15s solver is a good first choice for most stiff problems. However, the other stiff solvers might be more efficient for certain types of problems. This example solves a stiff test equation using all four stiff ODE solvers.
Consider the test equation
The equation becomes increasingly stiff as the magnitude of
and the initial condition
over the time interval
[0 0.5]. With these values, the problem is stiff enough that
ode23 struggle to integrate the equation. Also, use
odeset to pass in the constant Jacobian
and turn on the display of solver statistics.
lambda = 1e9; opts = odeset('Stats','on','Jacobian',-lambda); tspan = [0 0.5]; y0 = 1; subplot(2,2,1) disp('ode15s stats:') tic, ode15s(@(t,y) -lambda*y, tspan, y0, opts), toc title('ode15s') subplot(2,2,2) disp(' ') disp('ode23s stats:') tic, ode23s(@(t,y) -lambda*y, tspan, y0, opts), toc title('ode23s') subplot(2,2,3) disp(' ') disp('ode23t stats:') tic, ode23t(@(t,y) -lambda*y, tspan, y0, opts), toc title('ode23t') subplot(2,2,4) disp(' ') disp('ode23tb stats:') tic, ode23tb(@(t,y) -lambda*y, tspan, y0, opts), toc title('ode23tb')
ode15s stats: 104 successful steps 1 failed attempts 212 function evaluations 0 partial derivatives 21 LU decompositions 210 solutions of linear systems Elapsed time is 1.013618 seconds. ode23s stats: 63 successful steps 0 failed attempts 191 function evaluations 0 partial derivatives 63 LU decompositions 189 solutions of linear systems Elapsed time is 0.261956 seconds. ode23t stats: 95 successful steps 0 failed attempts 125 function evaluations 0 partial derivatives 28 LU decompositions 123 solutions of linear systems Elapsed time is 0.414460 seconds. ode23tb stats: 71 successful steps 0 failed attempts 167 function evaluations 0 partial derivatives 23 LU decompositions 236 solutions of linear systems Elapsed time is 0.308677 seconds.
The stiff solvers all perform well, but
ode23s completes the integration with the fewest steps and runs the fastest for this particular problem. Since the constant Jacobian is specified, none of the solvers need to calculate partial derivatives to compute the solution. Specifying the Jacobian benefits
ode23s the most since it normally evaluates the Jacobian in every step.
For general stiff problems, the performance of the stiff solvers varies depending on the format of the problem and specified options. Providing the Jacobian matrix or sparsity pattern always improves solver efficiency for stiff problems. But since the stiff solvers use the Jacobian differently, the improvement can vary significantly. Practically speaking, if a system of equations is very large or needs to be solved many times, then it is worthwhile to investigate the performance of the different solvers to minimize execution time.
The van der Pol equation is a second order ODE
Solve the van der Pol equation with
ode15s. The function
vdp1000.m ships with MATLAB® and encodes the equations. Specify a single output to return a structure containing information about the solution, such as the solver and evaluation points.
tspan = [0 3000]; y0 = [2 0]; sol = ode15s(@vdp1000,tspan,y0)
sol = solver: 'ode15s' extdata: [1x1 struct] x: [1x592 double] y: [2x592 double] stats: [1x1 struct] idata: [1x1 struct]
linspace to generate 2500 points in the interval
[0 3000]. Evaluate the first component of the solution at these points using
x = linspace(0,3000,2500); y = deval(sol,x,1);
Plot the solution.
Extend the solution to
odextend and add the result to the original plot.
tf = 4000; sol_new = odextend(sol,@vdp1000,tf); x = linspace(3000,tf,350); y = deval(sol_new,x,1); hold on plot(x,y,'r')
This example reformulates a system of ODEs as a system of differential algebraic equations (DAEs). The Robertson problem found in hb1ode.m is a classic test problem for programs that solve stiff ODEs. The system of equations is
hb1ode solves this system of ODEs to steady state with the initial conditions
. But the equations also satisfy a linear conservation law,
In terms of the solution and initial conditions, the conservation law is
The system of equations can be rewritten as a system of DAEs by using the conservation law to determine the state of . This reformulates the problem as the DAE system
The differential index of this system is 1, since only a single derivative of is required to make this a system of ODEs. Therefore, no further transformations are required before solving the system.
robertsdae encodes this DAE system. Save
robertsdae.m in your current folder to run the example.
function out = robertsdae(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 ];
The full example code for this formulation of the Robertson problem is available in hb1dae.m.
Solve the DAE system using
ode15s. Consistent initial conditions for
y0 are obvious based on the conservation law. Use
odeset to set the options:
Use a constant mass matrix to represent the left hand side of the system of equations.
Set the relative error tolerance to
Use an absolute tolerance of
1e-10 for the second solution component, since the scale varies dramatically from the other components.
'MassSingular' option at its default value
'maybe' to test the automatic detection of a DAE.
y0 = [1; 0; 0]; tspan = [0 4*logspace(-6,6)]; M = [1 0 0; 0 1 0; 0 0 0]; options = odeset('Mass',M,'RelTol',1e-4,'AbsTol',[1e-6 1e-10 1e-6]); [t,y] = ode15s(@robertsdae,tspan,y0,options);
Plot the solution.
y(:,2) = 1e4*y(:,2); semilogx(t,y); ylabel('1e4 * y(:,2)'); title('Robertson DAE problem with a Conservation Law, solved by ODE15S');
odefun— Functions to solvefunction handle
Functions to solve, specified as a function handle which defines the functions to be integrated.
dydt = odefun(t,y), for a scalar
a column vector
y, must return a column vector
corresponds to .
accept both input arguments,
even if one of the arguments is not used in the function.
For example, to solve , use the function:
function dydt = odefun(t,y) dydt = 5*y-3;
For a system of equations, the output of
a vector. Each element in the vector is the solution to one equation.
For example, to solve
use the function:
function dydt = odefun(t,y) dydt = zeros(2,1); dydt(1) = y(1)+2*y(2); dydt(2) = 3*y(1)+2*y(2);
For information on how to provide additional parameters to the
odefun, see Parameterizing Functions.
Interval of integration, specified as a vector. At minimum,
be a two element vector
[t0 tf] specifying the
initial and final times. To obtain solutions at specific times between
use a longer vector of the form
The elements in
tspan must be all increasing or
The solver imposes the initial conditions,
tspan(1), then integrates from
tspan has two elements,
tf], then the solver returns the solution evaluated at each
internal integration step within the interval.
tspan contains more than two
[t0,t1,t2,...,tf], then the solver returns
the solution evaluated at the given points. This does not affect
the internal steps that the solver uses to traverse from
Thus, the solver does not necessarily step precisely to each point
tspan. However, the solutions produced
at the specified points are of the same order of accuracy as the solutions
computed at each internal step.
Specifying several intermediate points has little effect on the efficiency of computation, but for large systems it can affect memory management.
The solution obtained by the solver might be different depending
on whether you specify
tspan as a two-element vector
or as a vector with intermediate points. If
several intermediate points, then they give an indication of the scale
for the problem, which can affect the size of the initial step taken
by the solver.
3 5 7 9 10]
y0— Initial conditionsvector
Initial conditions, specified as a vector.
be the same length as the vector output of
y0 contains an initial condition for each
equation defined in
options— Option structurestructure array
options = odeset('RelTol',1e-5,'Stats','on','OutputFcn',@odeplot) specifies
a relative error tolerance of
1e-5, turns on the
display of solver statistics, and specifies the output function
plot the solution as it is computed.
t— Evaluation pointscolumn vector
Evaluation points, returned as a column vector.
tspan contains two elements,
t contains the internal evaluation
points used to perform the integration.
tspan contains more than two
t is the same as
Solutions, returned as an array. Each row in
to the solution at the value returned in the corresponding row of
te— Time of eventscolumn vector
Time of events, returned as a column vector. The event times
te correspond to the solutions returned in
ie specifies which event occurred.
Solution at time of events, returned as an array. The event
te correspond to the solutions returned
ie specifies which
ie— Index of vanishing event functioncolumn vector
Index of vanishing event function, returned as a column vector.
The event times in
te correspond to the solutions
which event occurred.
sol— Structure for evaluationstructure array
Structure for evaluation, returned as a structure array. Use
this structure with the
to evaluate the solution at any point in the interval
sol structure array always includes
Row vector of the steps chosen by the solver.
Solutions. Each column
Additionally, if you specify
Events option and events are detected, then
includes these fields:
Points when events occurred.
Solutions that correspond to events in
Indices into the vector returned by the function specified
ode15s is a variable-step, variable-order
(VSVO) solver based on the numerical differentiation formulas (NDFs)
of orders 1 to 5. Optionally, it can use the backward differentiation
formulas (BDFs, also known as Gear's method) that are usually less
a multistep solver. Use
or is very inefficient and you suspect that the problem is stiff,
or when solving a differential-algebraic equation (DAE) , .