Differential Equations (Mathematics)

Mathematics

Changing ODE Integration Properties

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 one or more property values in an options structure.

[t,y] = solver(odefun,tspan,y0,options)

This section:

Explains how to create, modify, and query an options structure.
Describes the properties that you can use in an options structure.

In this and subsequent property tables, the most commonly used properties are listed first, followed by more advanced properties.

ODE Property Categories
Properties Category
Property Name

Error control
RelTol, AbsTol, NormControl

Solver output
OutputFcn, OutputSel, Refine, Stats

Jacobian matrix
Jacobian, JPattern, Vectorized

Step-size
InitialStep, MaxStep

Mass matrix and DAEs
Mass, MStateDependence, MvPattern, MassSingular, InitialSlope

Event location
Events

ode15s-specific
MaxOrder, BDF

**ODE Property Categories**
Properties Category	Property Name
Error control	`RelTol`, `AbsTol`, `NormControl`
Solver output	`OutputFcn`, `OutputSel`, `Refine`, `Stats`
Jacobian matrix	`Jacobian`, `JPattern`, `Vectorized`
Step-size	`InitialStep`, `MaxStep`
Mass matrix and DAEs	`Mass`, `MStateDependence`, `MvPattern`, `MassSingular`, `InitialSlope`
Event location	`Events`
`ode15s`-specific	`MaxOrder`, `BDF`

Creating and Maintaining an ODE Options Structure

Creating an Options Structure. The odeset function creates an options structure that you can pass as an argument to any of the ODE solvers. To create an options structure, odeset accepts property name/property value pairs using the syntax

options = odeset('name1',value1,'name2',value2,...)

In the resulting structure, options, the named properties have the specified values. Any unspecified properties contain default values in the solvers. For all properties, it is sufficient to type only the leading characters that uniquely identify the property name. odeset ignores case for property names.

With no input arguments, odeset displays all property names and their possible values. It indicates defaults with {}.

Modifying an Existing Options Structure. To modify an existing options structure, oldopts, use

options = odeset(oldopts,'name1',value1,...)

This sets options equal to the existing structure oldopts, overwrites any values in oldopts that are respecified using name/value pairs, and adds any new pairs to the structure. The modified structure is returned as an output argument. In the same way, the command

```
options = odeset(oldopts,newopts)
```

combines the structures oldopts and newopts. In the output argument, any values in the second argument overwrite those in the first argument.

Querying Options. The odeget function extracts property values from an options structure created with odeset.

```
o = odeget(options,'name')
```

This functions returns the value of the specified property, or an empty matrix [], if the property value is unspecified in the options structure.

As with odeset, it is sufficient to type only the leading characters that uniquely identify the property name. Case is ignored for property names.

Error Control Properties

At each step, the solver estimates the local error e in the ith component of the solution. This error must be less than or equal to the acceptable error, which is a function of the specified relative tolerance, RelTol, and the specified absolute tolerance, AbsTol.

|e(i)|  max(RelTol*abs(y(i)),AbsTol(i))

For routine problems, the ODE solvers deliver accuracy roughly equivalent to the accuracy you request. They deliver less accuracy for problems integrated over "long" intervals and problems that are moderately unstable. Difficult problems may require tighter tolerances than the default values. For relative accuracy, adjust RelTol. For the absolute error tolerance, the scaling of the solution components is important: if |y| is somewhat smaller than AbsTol, the solver is not constrained to obtain any correct digits in y. You might have to solve a problem more than once to discover the scale of solution components.

Roughly speaking, this means that you want RelTol correct digits in all solution components except those smaller than thresholds AbsTol(i). Even if you are not interested in a component y(i) when it is small, you may have to specify AbsTol(i) small enough to get some correct digits in y(i) so that you can accurately compute more interesting components

The following table describes the error control properties. Use odeset to set the properties.

ODE Error Control Properties
Property
Value
Description

RelTol
Positive scalar {1e-3}
A relative error tolerance that applies to all components of the solution vector y. It is a measure of the error relative to the size of each solution component. Roughly, it controls the number of correct digits in all solution components except those smaller than thresholds AbsTol(i).
The default, 1e-3, corresponds to 0.1% accuracy.

AbsTol
Positive scalar or vector {1e-6}
Absolute error tolerances that apply to the individual components of the solution vector. AbsTol(i) is a threshold below which the value of the ith solution component is unimportant. The absolute error tolerances determine the accuracy when the solution approaches zero. Even if you are not interested in a component y(i) when it is small, you may have to specify AbsTol(i) small enough to get some correct digits in y(i) so that you can accurately compute more interesting components.
If AbsTol is a vector, the length of AbsTol must be the same as the length of the solution vector y. If AbsTol is a scalar, the value applies to all components of y.

NormControl
on | {off}
Control error relative to norm of solution. Set this property on to request that the solvers control the error in each integration step with norm(e) <= max(RelTol*norm(y),AbsTol). By default the solvers use a more stringent component-wise error control.

**ODE Error Control Properties**
Property	Value	Description
`RelTol`	Positive scalar {`1e-3`}	A relative error tolerance that applies to all components of the solution vector `y`. It is a measure of the error relative to the size of each solution component. Roughly, it controls the number of correct digits in all solution components except those smaller than thresholds `AbsTol(i)`. The default, `1e-3`, corresponds to 0.1% accuracy.
`AbsTol`	Positive scalar or vector {`1e-6`}	Absolute error tolerances that apply to the individual components of the solution vector. `AbsTol(i)` is a threshold below which the value of the `i`th solution component is unimportant. The absolute error tolerances determine the accuracy when the solution approaches zero. Even if you are not interested in a component `y(i)` when it is small, you may have to specify `AbsTol(i)` small enough to get some correct digits in `y(i)` so that you can accurately compute more interesting components. If `AbsTol` is a vector, the length of `AbsTol` must be the same as the length of the solution vector `y`. If `AbsTol` is a scalar, the value applies to all components of `y`.
`NormControl`	`on` \| {`off`}	Control error relative to norm of solution. Set this property `on` to request that the solvers control the error in each integration step with `norm(e) <= max(RelTol*norm(y),AbsTol)`. By default the solvers use a more stringent component-wise error control.

Solver Output Properties

The solver output properties let you control the output that the solvers generate. Use odeset to set these properties.

ODE Solver Output Properties
Property
Value
Description

OutputFcn
Function {odeplot}
Installable output function. The solver calls this function after every successful integration step.
For example,
options = odeset('OutputFcn',@myfun)
sets the OutputFcn property to an output function, myfun, that can be passed to an ODE solver.
The output function must be of the form
status = myfun(t,y,flag,p1,p2,...)
The solver calls the specified output function with the following flags. Note that the syntax of the call differs with the flag. The function must respond appropriately:

init
The solver calls myfun(tspan,y0,'init') before beginning the integration to allow the output function to initialize. tspan and y0 are the input arguments to the ODE solver.

{none}
The solver calls status = myfun(t,y) after each integration step on which output is requested. t contains points where output was generated during the step, and y is the numerical solution at the points in t. If t is a vector, the ith column of y corresponds to the ith element of t.
When length(tspan) > 2 the output is produced at every point in tspan. When length(tspan) = 2 the output is produced according to the Refine option.
myfun must return a status output value of 0 or 1. If status = 1, the solver halts integration. You can use this mechanism, for instance, to implement a Stop button.

done
The solver calls myfun([],[],'done') when integration is complete to allow the output function to perform any cleanup chores.

You can use these general purpose output functions or you can edit them to create your own. Type help function at the command line for more information.
odeplot - time series plotting (default when you call the solver with no output arguments and you have not specified an output function)

odephas2 - two-dimensional phase plane plotting

odephas3 - three-dimensional phase plane plotting

odeprint - print solution as it is computed

Note If you call the solver with no output arguments, the solver does not allocate storage to hold the entire solution history.

OutputSel
Vector of indices
Vector of indices specifying which components of the solution vector are to be passed to the output function. For example, if you want to use the odeplot output function, but you want to plot only the first and third components of the solution, you can do this using
options = odeset('OutputFcn',@odeplot,'OutputSel',[1 3]);

By default, the solver passes all components of the solution to the output function.

Refine
Positive integer
Increases the number of output points by a factor of Refine. If Refine is 1, the solver returns solutions only at the end of each time step. If Refine is n >1, the solver subdivides each time step into n smaller intervals, and returns solutions at each time point. Refine does not apply when length(tspan)>2.
Note In all the solvers, the default value of Refine is 1. Within ode45, however, the default is 4 to compensate for the solver's large step sizes. To override this and see only the time steps chosen by ode45, set Refine to 1.
The extra values produced for Refine are computed by means of continuous extension formulas. These are specialized formulas used by the ODE solvers to obtain accurate solutions between computed time steps without significant increase in computation time.

Stats
on | {off}
Specifies whether the solver should display statistics about its computations. By default, Stats is off. If it is on, after solving the problem the solver displays:
The number of successful steps

The number of failed attempts

The number of times the ODE function was called to evaluate

The number of times that the partial derivatives matrix was formed

The number of LU decompositions

The number of solutions of linear systems

Jacobian Matrix Properties

The stiff ODE solvers often execute faster if you provide additional information about the Jacobian matrix , a matrix of partial derivatives of the function that defines the differential equations.

The Jacobian matrix properties pertain only to those solvers for stiff problems (ode15s, ode23s, ode23t, and ode23tb) for which the Jacobian matrix can be critical to reliability and efficiency. If you do not provide a function to calculate the Jacobian, these solvers approximate the Jacobian numerically using finite differences. In this case, you may want to use the Vectorized, or JPattern properties.

The following table describes the Jacobian matrix properties. Use odeset to set these properties.

ODE Jacobian Matrix Properties
Property
Value
Description

Jacobian
Function | constant matrix
A constant matrix or a function that evaluates the Jacobian. Supplying an analytical Jacobian often increases the speed and reliability of the solution for stiff problems. Set this property to a function FJac, where FJac(t,y) computes , or to the constant value of .
The Jacobian for the stiff van der Pol problem shown above can be coded as
function J = vdp1000jac(t,y) J = [ 0 1 (-2000*y(1)*y(2)-1) (1000*(1-y(1)^2)) ];

JPattern
Sparse matrix of {0,1}
Sparsity pattern with 1s where there might be nonzero entries in the Jacobian. It is used to generate a sparse Jacobian matrix numerically.
Set this property to a sparse matrix with if component of depends on component of , and 0 otherwise. The solver uses this sparsity pattern to generate a sparse Jacobian matrix numerically. If the Jacobian matrix is large and sparse, this can greatly accelerate execution. For an example using the JPattern property, see Example: Large, Stiff Sparse Problem (brussode).

Vectorized
on | {off}
Set on to inform the solver that you have coded the ODE function F so that F(t,[y1 y2 ...]) returns [F(t,y1) F(t,y2) ...]. This allows the solver to reduce the number of function evaluations required to compute all the columns of the Jacobian matrix, and may significantly reduce solution time.
With the MATLAB array notation, it is typically an easy matter to vectorize an ODE function. For example, the stiff van der Pol example shown previously can be vectorized by introducing colon notation into the subscripts and by using the array power (.^) and array multiplication (.*) operators.
function dydt = vdp1000(t,y) dydt = [y(2,:); 1000*(1-y(1,:).^2).*y(2,:)-y(1,:)];

Note Vectorization of the ODE function used by the ODE solvers differs from the vectorization used by the BVP solver, bvp4c. For the ODE solvers, the ODE function is vectorized only with respect to the second argument, while bvp4c requires vectorization with respect the first and second arguments.

Step-Size Properties

The step-size properties let you specify the size of the first step the solver tries, potentially helping it to better recognize the scale of the problem. In addition, you can specify bounds on the sizes of subsequent time steps.

The following table describes the step-size properties. Use odeset to set these properties.

ODE Step Size Properties
Property
Value
Description

InitialStep
Positive scalar
Suggested initial step size. InitialStep sets an upper bound on the magnitude of the first step size the solver tries. If you do not set InitialStep, the initial step size is based on the slope of the solution at the initial time tspan(1), and if the slope of all solution components is zero, the procedure might try a step size that is much too large. If you know this is happening or you want to be sure that the solver resolves important behavior at the start of the integration, help the code start by providing a suitable InitialStep.

MaxStep
Positive scalar {0.1*abs(t0-tf)}
Upper bound on solver step size. If the differential equation has periodic coefficients or solutions, it may be a good idea to set MaxStep to some fraction (such as 1/4) of the period. This guarantees that the solver does not enlarge the time step too much and step over a period of interest. Do not reduce MaxStep:
To produce more output points. This can significantly slow down solution time. Instead, use Refine to compute additional outputs by continuous extension at very low cost.

When the solution does not appear to be accurate enough. Instead, reduce the relative error tolerance RelTol, and use the solution you just computed to determine appropriate values for the absolute error tolerance vector AbsTol. (See Error Control Properties for a description of the error tolerance properties.)

To make sure that the solver doesn't step over some behavior that occurs only once during the simulation interval. If you know the time at which the change occurs, break the simulation interval into two pieces and call the solvers twice. If you do not know the time at which the change occurs, try reducing the error tolerances RelTol and AbsTol. Use MaxStep as a last resort.

**ODE Step Size Properties**
Property	Value	Description
`InitialStep`	Positive scalar	Suggested initial step size. `InitialStep` sets an upper bound on the magnitude of the first step size the solver tries. If you do not set `InitialStep`, the initial step size is based on the slope of the solution at the initial time `tspan(1)`, and if the slope of all solution components is zero, the procedure might try a step size that is much too large. If you know this is happening or you want to be sure that the solver resolves important behavior at the start of the integration, help the code start by providing a suitable `InitialStep`.
`MaxStep`	Positive scalar {`0.1`*`abs(t0-tf)`}	Upper bound on solver step size. If the differential equation has periodic coefficients or solutions, it may be a good idea to set `MaxStep` to some fraction (such as 1/4) of the period. This guarantees that the solver does not enlarge the time step too much and step over a period of interest. Do not reduce `MaxStep`: To produce more output points. This can significantly slow down solution time. Instead, use `Refine` to compute additional outputs by continuous extension at very low cost. When the solution does not appear to be accurate enough. Instead, reduce the relative error tolerance `RelTol`, and use the solution you just computed to determine appropriate values for the absolute error tolerance vector `AbsTol`. (See Error Control Properties for a description of the error tolerance properties.)
		To make sure that the solver doesn't step over some behavior that occurs only once during the simulation interval. If you know the time at which the change occurs, break the simulation interval into two pieces and call the solvers twice. If you do not know the time at which the change occurs, try reducing the error tolerances `RelTol` and `AbsTol`. Use `MaxStep` as a last resort.

Mass Matrix and DAE Properties

The solvers of the ODE suite can solve ODEs of the form

(14-2)

with a mass matrix that can be sparse.

When is nonsingular, the equation above is equivalent to and the ODE has a solution for any initial values at . The more general form (Equation 14-2) is convenient when you express a model naturally in terms of a mass matrix. For large, sparse , solving Equation 14-2 directly reduces the storage and runtime needed to solve the problem.

When is singular, then is a differential-algebraic equation (DAE). A DAE has a solution only when is consistent, that is, there exists an initial slope such that

. If and 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. For DAEs of index 1, solving an initial value problem with consistent initial conditions is much like solving an ODE.

The ode15s and ode23t solvers can solve DAEs of index 1. For examples of DAE problems, see hb1dae (Example: Differential-Algebraic Problem) and amp1dae.

The following table describes the mass matrix and DAE properties. Use odeset to set these properties.

ODE Mass Matrix and DAE Properties
Property
Value
Description

Mass
Constant matrix | function
Constant mass matrix or a function that evaluates the mass matrix . For problems set this property to the value of the constant mass matrix . For problems , set this property to a function Mfun, where Mfun(t,y) evaluates the mass matrix . When solving DAEs, it is advantageous to formulate the problem so that is a diagonal matrix (a semi-explicit DAE). The ode23s solver can only solve problems with a constant mass matrix .
For example problems, see fem1ode (Example: Finite Element Discretization), fem2ode, or batonode.

MStateDependence
none | {weak} | strong
Dependence of the mass matrix on . Set this property to none for problems . Both weak and strong indicate , but weak results in implicit solvers using approximations when solving algebraic equations.

MvPattern
Sparse matrix
sparsity pattern. Set this property to a sparse matrix with if for any , the component of depends on component of , and 0 otherwise. For use with the ode15s, ode23t, and ode23tb solvers when MStateDependence is strong. See burgersode as an example.

MassSingular
yes | no | {maybe}
Indicates whether the mass matrix is singular. Set this property to no if the mass matrix is not singular and you are using either the ode15s or ode23t solver. The default value of maybe causes the solver to test whether the problem is a DAE, i.e., whether is singular.
For an example of a problem with a mass matrix, see Example: Finite Element Discretization (fem1ode).

InitialSlope
Vector | {zero vector}
Vector representing the consistent initial slope , where satisfies

. The default is the zero vector.
This property is for use with the ode15s and ode23t solvers when solving DAEs.

**ODE Mass Matrix and DAE Properties**
Property	Value	Description
`Mass`	Constant matrix \| function	Constant mass matrix or a function that evaluates the mass matrix . For problems set this property to the value of the constant mass matrix . For problems , set this property to a function `Mfun`, where `Mfun(t,y)` evaluates the mass matrix . When solving DAEs, it is advantageous to formulate the problem so that is a diagonal matrix (a semi-explicit DAE). The `ode23s` solver can only solve problems with a constant mass matrix . For example problems, see `fem1ode` (Example: Finite Element Discretization), `fem2ode`, or `batonode`.
`MStateDependence`	`none \| {weak} \| strong`	Dependence of the mass matrix on . Set this property to `none` for problems . Both `weak` and `strong` indicate , but `weak` results in implicit solvers using approximations when solving algebraic equations.
`MvPattern`	Sparse matrix	sparsity pattern. Set this property to a sparse matrix with if for any , the component of depends on component of , and 0 otherwise. For use with the `ode15s`, `ode23t`, and `ode23tb` solvers when `MStateDependence` is `strong`. See `burgersode` as an example.
`MassSingular`	`yes \| no \| {maybe}`	Indicates whether the mass matrix is singular. Set this property to `no` if the mass matrix is not singular and you are using either the `ode15s` or `ode23t` solver. The default value of `maybe` causes the solver to test whether the problem is a DAE, i.e., whether is singular. For an example of a problem with a mass matrix, see Example: Finite Element Discretization (`fem1ode`).
`InitialSlope`	Vector \| {zero vector}	Vector representing the consistent initial slope , where satisfies . The default is the zero vector. This property is for use with the `ode15s` and `ode23t` solvers when solving DAEs.

Event Location Property

In some ODE problems the times of specific events are important, such as the time at which a ball hits the ground, or the time at which a spaceship returns to the earth. While solving a problem, the ODE solvers can detect such events by locating transitions to, from, or through zeros of user-defined functions.

The following table describes the Events property. Use odeset to set this property.

ODE Events Property
String
Value
Description

Events
Function
Function that includes one or more event functions. The function is of the form
[value,isterminal,direction] = events(t,y)
value, isterminal, and direction are vectors for which the ith element corresponds to the ith event function:
value(i) is the value of the ith event 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 located (the default), +1 if only zeros where the event function is increasing, and -1 if only zeros where the event function is decreasing.

If you specify an events function and events are detected, the solver returns three additional outputs:
A column vector of times at which events occur

Solution values corresponding to these times

Indices into the vector returned by the events function. The values indicate which event the solver detected.

If you call the solver as
[T,Y,TE,YE,IE] = solver(odefun,tspan,y0,options)
the solver returns these outputs as TE, YE, and IE respectively. If you call the solver as
sol = solver(odefun,tspan,y0,options)

the solver returns these outputs as sol.xe, sol.ye, and sol.ie respectively.

For examples that use an event function, see Example: Simple Event Location (ballode) and Example: Advanced Event Location (orbitode).

**ODE Events Property**
String	Value	Description
`Events`	Function	Function that includes one or more event functions. The function is of the form [`value`,`isterminal`,`direction`] = `events`(`t`,`y`) `value`, `isterminal`, and `direction` are vectors for which the `i`th element corresponds to the `i`th event function: `value(i)` is the value of the `i`th event 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 located (the default), `+1` if only zeros where the event function is increasing, and `-1` if only zeros where the event function is decreasing. If you specify an events function and events are detected, the solver returns three additional outputs: `A` column vector of times at which events occur Solution values corresponding to these times Indices into the vector returned by the events function. The values indicate which event the solver detected. If you call the solver as [T,Y,TE,YE,IE] = `solver`(odefun,tspan,y0,options) the solver returns these outputs as `TE`, `YE`, and `IE` respectively. If you call the solver as sol = `solver`(odefun,tspan,y0,options) the solver returns these outputs as `sol.xe`, `sol.ye`, and `sol.ie` respectively. For examples that use an event function, see Example: Simple Event Location (`ballode`) and Example: Advanced Event Location (`orbitode`).

ode15s Properties

ode15s is a variable-order solver for stiff problems. It is based on the numerical differentiation formulas (NDFs). The NDFs are generally more efficient than the closely related family of backward differentiation formulas (BDFs), also known as Gear's methods. The ode15s properties let you choose between these formulas, as well as specifying the maximum order for the formula used.

The following table describes the ode15s properties. Use odeset to set these properties.

ode15s Properties
Property
Value
Description

MaxOrder
1 | 2 | 3 | 4 | {5}
The maximum order formula used to compute the solution.

BDF
on | {off}
Specifies whether you want to use the BDFs instead of the default NDFs. Set BDF on to have ode15s use the BDFs.
For both the NDFs and BDFs, the formulas of orders 1 and 2 are A-stable (the stability region includes the entire left half complex plane). The higher order formulas are not as stable, and the higher the order the worse the stability. There is a class of stiff problems (stiff oscillatory) that is solved more efficiently if MaxOrder is reduced (for example to 2) so that only the most stable formulas are used.

**ode15s Properties**
Property	Value	Description
`MaxOrder`	`1` \| `2` \| `3` \| `4` \| `{5}`	The maximum order formula used to compute the solution.
`BDF`	`on` \| `{off}`	Specifies whether you want to use the BDFs instead of the default NDFs. Set `BDF` `on` to have `ode15s` use the BDFs. For both the NDFs and BDFs, the formulas of orders 1 and 2 are A-stable (the stability region includes the entire left half complex plane). The higher order formulas are not as stable, and the higher the order the worse the stability. There is a class of stiff problems (stiff oscillatory) that is solved more efficiently if `MaxOrder` is reduced (for example to `2`) so that only the most stable formulas are used.

Solving ODE Problems Examples: Applying the ODE Initial Value Problem Solvers