Solve stiff differential equations — trapezoidal rule + backward differentiation formula
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, .
The solvers all use similar syntaxes. The
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.
a structure that you can use with
sol = ode23tb(___)
deval to evaluate
the solution at any point on the interval
You can use any of the input argument combinations in previous syntaxes.
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] = ode23tb(@(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 and . 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
ode23tb solver, and then plot the first column of the solution
y against the time points
ode23tb solver passes through stiff areas with far fewer steps than
[t,y] = ode23tb(@vdp1000,[0 3000],[2 0]); plot(t,y(:,1),'-o')
ode23s 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
ode23tb. 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] = ode23tb(@(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 increases. Use 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; y0 = 1; tspan = [0 0.5]; opts = odeset('Jacobian',-lambda,'Stats','on');
Solve the equation with
ode23tb. Make subplots for comparison.
subplot(2,2,1) tic, ode15s(@(t,y) -lambda*y, tspan, y0, opts), toc
104 successful steps 1 failed attempts 212 function evaluations 0 partial derivatives 21 LU decompositions 210 solutions of linear systems Elapsed time is 1.479842 seconds.
title('ode15s') subplot(2,2,2) tic, ode23s(@(t,y) -lambda*y, tspan, y0, opts), toc
63 successful steps 0 failed attempts 191 function evaluations 0 partial derivatives 63 LU decompositions 189 solutions of linear systems Elapsed time is 0.408763 seconds.
title('ode23s') subplot(2,2,3) tic, ode23t(@(t,y) -lambda*y, tspan, y0, opts), toc
95 successful steps 0 failed attempts 125 function evaluations 0 partial derivatives 28 LU decompositions 123 solutions of linear systems Elapsed time is 0.674248 seconds.
title('ode23t') subplot(2,2,4) tic, ode23tb(@(t,y) -lambda*y, tspan, y0, opts), toc
71 successful steps 0 failed attempts 167 function evaluations 0 partial derivatives 23 LU decompositions 236 solutions of linear systems Elapsed time is 0.562609 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.
ode23tb is an implementation of TR-BDF2,
an implicit Runge-Kutta formula with a trapezoidal rule step as its
first stage and a backward differentiation formula of order two as
its second stage. By construction, the same iteration matrix is used
in evaluating both stages. Like
this solver may be more efficient than
problems with crude tolerances , .
 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.
 Shampine, L. F. and M. E. Hosea, “Analysis and Implementation of TR-BDF2,” Applied Numerical Mathematics 20, 1996.