Documentation

This is machine translation

Translated by Microsoft
Mouse over text to see original. Click the button below to return to the English verison of the page.

ode23s

Solve stiff differential equations — low order method

Syntax

  • [t,y] = ode23s(odefun,tspan,y0)
    example
  • [t,y] = ode23s(odefun,tspan,y0,options)
    example
  • [t,y,te,ye,ie] = ode23s(odefun,tspan,y0,options)
  • sol = ode23s(___)

Description

example

[t,y] = ode23s(odefun,tspan,y0), where tspan = [t0 tf], integrates the system of differential equations y'=f(t,y) from t0 to tf with initial conditions y0. Each row in the solution array y corresponds to a value returned in column vector t.

All MATLAB® ODE solvers can solve systems of equations of the form y'=f(t,y), or problems that involve a mass matrix, M(t,y)y'=f(t,y). The ode23s solver only can solve problems with a mass matrix if the mass matrix is constant. ode15s and ode23t can solve problems with a mass matrix that is singular, known as differential-algebraic equations (DAEs). Specify the mass matrix using the Mass option of odeset.

example

[t,y] = ode23s(odefun,tspan,y0,options) also uses the integration settings defined by options, which is an argument created using the odeset function. For example, use the AbsTol and RelTol options to specify absolute and relative error tolerances, or the Mass option to provide a mass matrix.

[t,y,te,ye,ie] = ode23s(odefun,tspan,y0,options) additionally finds where functions of (t,y), called event functions, are zero. In the output, te is the time of the event, ye is the solution at the time of the event, and ie is the index of the triggered event.

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 'Events' property to a function, such as myEventFcn or @myEventFcn, and creating a corresponding function: [value,isterminal,direction] = myEventFcn(t,y). For more information, see ODE Event Location.

sol = ode23s(___) returns a structure that you can use with deval to evaluate the solution at any point on the interval [t0 tf]. You can use any of the input argument combinations in previous syntaxes.

Examples

collapse all

ODE With Single Solution Component

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

$$y' = -10t.$$

Use a time interval of [0,2] and the initial condition y0 = 1.

tspan = [0 2];
y0 = 1;
[t,y] = ode23s(@(t,y) -10*t, tspan, y0);

Plot the solution.

plot(t,y,'-o')

Solve Stiff ODE

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:

$$\begin{array}{cl} y_1' &= y_2\\y_2' &= 1000(1-y_1^2)y_2-y_1\end{array}$$

The initial conditions are $y_1(0)=2$ and $y_2(0)=0$. 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-3 and 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 ode23s solver, and then plot the first column of the solution y against the time points t. The ode23s solver passes through stiff areas with far fewer steps than ode45.

[t,y] = ode23s(@vdp1000,[0 3000],[2 0]);
plot(t,y(:,1),'-o')

Pass Extra Parameters to ODE Function

ode23s only works with functions that use two input arguments, t and 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

$$y'' = \frac{A}{B} t y.$$

Rewriting the equation as a first-order system yields

$$\begin{array}{cl} y'_1 &= y_2\\ y'_2 &= \frac{A}{B} t y_1.
\end{array}$$

odefcn.m represents this system of equations as a function that accepts four input arguments: t, y, A, and B.

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 ode23s. Specify the function handle such that it passes in the predefined values for A and B to odefcn.

A = 1;
B = 2;
tspan = [0 5];
y0 = [0 0.01];
[t,y] = ode23s(@(t,y) odefcn(t,y,A,B), tspan, y0);

Plot the results.

plot(t,y(:,1),'-o',t,y(:,2),'-.')

Compare Stiff ODE Solvers

The 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

$$y' = -\lambda y.$$

The equation becomes increasingly stiff as the magnitude of $\lambda$ increases. Use $\lambda = 1 \times 10^9$ and the initial condition $y(0)=1$ over the time interval [0 0.5]. With these values, the problem is stiff enough that ode45 and ode23 struggle to integrate the equation. Also, use odeset to pass in the constant Jacobian $J = \frac{\partial f}{\partial y} = - \lambda$ 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.

Input Arguments

collapse all

odefun — Functions to solvefunction handle

Functions to solve, specified as a function handle which defines the functions to be integrated.

The function dydt = odefun(t,y), for a scalar t and a column vector y, must return a column vector dydt of data type single or double that corresponds to f(t,y). odefun must accept both input arguments, t and y, even if one of the arguments is not used in the function.

For example, to solve y'=5y3, use the function:

function dydt = odefun(t,y)
dydt = 5*y-3;

For a system of equations, the output of odefun is a vector. Each element in the vector is the solution to one equation. For example, to solve

y'1=y1+2y2y'2=3y1+2y2

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 function odefun, see Parameterizing Functions.

Example: @myFcn

Data Types: function_handle

tspan — Interval of integrationvector

Interval of integration, specified as a vector. At minimum, tspan must be a two element vector [t0 tf] specifying the initial and final times. To obtain solutions at specific times between t0 and tf, use a longer vector of the form [t0,t1,t2,...,tf]. The elements in tspan must be all increasing or all decreasing.

The solver imposes the initial conditions, y0, at tspan(1), then integrates from tspan(1) to tspan(end):

  • If tspan has two elements, [t0 tf], then the solver returns the solution evaluated at each internal integration step within the interval.

  • If tspan contains more than two elements [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 tspan(1) to tspan(end). Thus, the solver does not necessarily step precisely to each point specified in 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 tspan contains 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.

Example: [1 10]

Example: [1 3 5 7 9 10]

Data Types: single | double

y0 — Initial conditionsvector

Initial conditions, specified as a vector. y0 must be the same length as the vector output of odefun, so that y0 contains an initial condition for each equation defined in odefun.

Data Types: single | double

options — Option structurestructure array

Option structure, specified as a structure array. Use the odeset function to create or modify the options structure. See Summary of ODE Options for a list of the options compatible with each solver.

Example: 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 @odeplot to plot the solution as it is computed.

Data Types: struct

Output Arguments

collapse all

t — Evaluation pointscolumn vector

Evaluation points, returned as a column vector.

  • If tspan contains two elements, [t0 tf], then t contains the internal evaluation points used to perform the integration.

  • If tspan contains more than two elements, then t is the same as tspan.

y — Solutionsarray

Solutions, returned as an array. Each row in y corresponds to the solution at the value returned in the corresponding row of t.

te — Time of eventscolumn vector

Time of events, returned as a column vector. The event times in te correspond to the solutions returned in ye, and ie specifies which event occurred.

ye — Solution at time of eventsarray

Solution at time of events, returned as an array. The event times in te correspond to the solutions returned in ye, and ie specifies which event occurred.

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 returned in ye, and ie specifies which event occurred.

sol — Structure for evaluationstructure array

Structure for evaluation, returned as a structure array. Use this structure with the deval function to evaluate the solution at any point in the interval [t0 tf]. The sol structure array always includes these fields:

Structure FieldDescription

sol.x

Row vector of the steps chosen by the solver.

sol.y

Solutions. Each column sol.y(:,i) contains the solution at time sol.x(i).

sol.solver

Solver name.

Additionally, if you specify the Events option and events are detected, then sol also includes these fields:

Structure FieldDescription

sol.xe

Points when events occurred. sol.xe(end) contains the exact point of a terminal event, if any.

sol.ye

Solutions that correspond to events in sol.xe.

sol.ie

Indices into the vector returned by the function specified in the Events option. The values indicate which event the solver detected.

More About

collapse all

Algorithms

ode23s is based on a modified Rosenbrock formula of order 2. Because it is a single-step solver, it may be more efficient than ode15s at solving problems that permit crude tolerances or problems with solutions that change rapidly. It can solve some kinds of stiff problems for which ode15s is not effective. The ode23s solver evaluates the Jacobian during each step of the integration, so supplying it with the Jacobian matrix is critical to its reliability and efficiency [1].

References

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

See Also

| | |

Introduced before R2006a

Was this topic helpful?