Previous: Differential-Algebraic Equations, Up: Differential Equations [Contents][Index]

Octave also provides a set of solvers for initial value problems for ordinary differential equations (ODEs) that have a MATLAB-compatible interface. The options for this class of methods are set using the functions.

Currently implemented solvers are:

- Runge-Kutta methods
- ode45 integrates a system of non-stiff ODEs or
index-1 differential-algebraic equations (DAEs) using the high-order,
variable-step Dormand-Prince method. It requires six function
evaluations per integration step, but may take larger steps on smooth
problems than
`ode23`

: potentially offering improved efficiency at smaller tolerances. - ode23 integrates a system of non-stiff ODEs or (or index-1 DAEs). It uses the third-order Bogacki-Shampine method and adapts the local step size in order to satisfy a user-specified tolerance. The solver requires three function evaluations per integration step.
- ode23s integrates a system of stiff ODEs (or index-1 DAEs) using a modified second-order Rosenbrock method.

- ode45 integrates a system of non-stiff ODEs or
index-1 differential-algebraic equations (DAEs) using the high-order,
variable-step Dormand-Prince method. It requires six function
evaluations per integration step, but may take larger steps on smooth
problems than
- Linear multistep methods
- ode15s integrates a system of stiff ODEs (or index-1 DAEs) using a variable step, variable order method based on Backward Difference Formulas (BDF).
- ode15i integrates a system of fully-implicit ODEs
(or index-1 DAEs) using the same variable step, variable order method as
`ode15s`

. decic can be used to compute consistent initial conditions for`ode15i`

.

Detailed information on the solvers are given in L. F. Shampine and M. W. Reichelt, The MATLAB ODE Suite, SIAM Journal on Scientific Computing, Vol. 18, 1997, pp. 1–22, DOI: https://doi.org/10.1137/S1064827594276424.

- :
*[*`t`,`y`] =**ode45***(*`fun`,`trange`,`init`) - :
*[*`t`,`y`] =**ode45***(*`fun`,`trange`,`init`,`ode_opt`) - :
*[*`t`,`y`,`te`,`ye`,`ie`] =**ode45***(…)* - :
`solution`=**ode45***(…)* - :
**ode45***(…)* -
Solve a set of non-stiff Ordinary Differential Equations (non-stiff ODEs) with the well known explicit Dormand-Prince method of order 4.

`fun`is a function handle, inline function, or string containing the name of the function that defines the ODE:`y' = f(t,y)`

. The function must accept two inputs where the first is time`t`and the second is a column vector of unknowns`y`.`trange`specifies the time interval over which the ODE will be evaluated. Typically, it is a two-element vector specifying the initial and final times (`[tinit, tfinal]`

). If there are more than two elements then the solution will also be evaluated at these intermediate time instances.By default,

`ode45`

uses an adaptive timestep with the`integrate_adaptive`

algorithm. The tolerance for the timestep computation may be changed by using the options`"RelTol"`

and`"AbsTol"`

.`init`contains the initial value for the unknowns. If it is a row vector then the solution`y`will be a matrix in which each column is the solution for the corresponding initial value in`init`.The optional fourth argument

`ode_opt`specifies non-default options to the ODE solver. It is a structure generated by`odeset`

.The function typically returns two outputs. Variable

`t`is a column vector and contains the times where the solution was found. The output`y`is a matrix in which each column refers to a different unknown of the problem and each row corresponds to a time in`t`.The output can also be returned as a structure

`solution`which has a field`x`containing a row vector of times where the solution was evaluated and a field`y`containing the solution matrix such that each column corresponds to a time in`x`. Use`fieldnames (`

to see the other fields and additional information returned.`solution`)If no output arguments are requested, and no

`"OutputFcn"`

is specified in`ode_opt`, then the`"OutputFcn"`

is set to`odeplot`

and the results of the solver are plotted immediately.If using the

`"Events"`

option then three additional outputs may be returned.`te`holds the time when an Event function returned a zero.`ye`holds the value of the solution at time`te`.`ie`contains an index indicating which Event function was triggered in the case of multiple Event functions.Example: Solve the Van der Pol equation

fvdp = @(

`t`,`y`) [`y`(2); (1 -`y`(1)^2) *`y`(2) -`y`(1)]; [`t`,`y`] = ode45 (fvdp, [0, 20], [2, 0]);

- :
*[*`t`,`y`] =**ode23***(*`fun`,`trange`,`init`) - :
*[*`t`,`y`] =**ode23***(*`fun`,`trange`,`init`,`ode_opt`) - :
*[*`t`,`y`,`te`,`ye`,`ie`] =**ode23***(…)* - :
`solution`=**ode23***(…)* - :
**ode23***(…)* -
Solve a set of non-stiff Ordinary Differential Equations (non-stiff ODEs) with the well known explicit Bogacki-Shampine method of order 3.

`fun`is a function handle, inline function, or string containing the name of the function that defines the ODE:`y' = f(t,y)`

. The function must accept two inputs where the first is time`t`and the second is a column vector of unknowns`y`.`trange`specifies the time interval over which the ODE will be evaluated. Typically, it is a two-element vector specifying the initial and final times (`[tinit, tfinal]`

). If there are more than two elements then the solution will also be evaluated at these intermediate time instances.By default,

`ode23`

uses an adaptive timestep with the`integrate_adaptive`

algorithm. The tolerance for the timestep computation may be changed by using the options`"RelTol"`

and`"AbsTol"`

.`init`contains the initial value for the unknowns. If it is a row vector then the solution`y`will be a matrix in which each column is the solution for the corresponding initial value in`init`.The optional fourth argument

`ode_opt`specifies non-default options to the ODE solver. It is a structure generated by`odeset`

.The function typically returns two outputs. Variable

`t`is a column vector and contains the times where the solution was found. The output`y`is a matrix in which each column refers to a different unknown of the problem and each row corresponds to a time in`t`.The output can also be returned as a structure

`solution`which has a field`x`containing a row vector of times where the solution was evaluated and a field`y`containing the solution matrix such that each column corresponds to a time in`x`. Use`fieldnames (`

to see the other fields and additional information returned.`solution`)If no output arguments are requested, and no

`"OutputFcn"`

is specified in`ode_opt`, then the`"OutputFcn"`

is set to`odeplot`

and the results of the solver are plotted immediately.If using the

`"Events"`

option then three additional outputs may be returned.`te`holds the time when an Event function returned a zero.`ye`holds the value of the solution at time`te`.`ie`contains an index indicating which Event function was triggered in the case of multiple Event functions.Example: Solve the Van der Pol equation

fvdp = @(

`t`,`y`) [`y`(2); (1 -`y`(1)^2) *`y`(2) -`y`(1)]; [`t`,`y`] = ode23 (fvdp, [0, 20], [2, 0]);Reference: For the definition of this method see https://en.wikipedia.org/wiki/List_of_Runge%E2%80%93Kutta_methods.

- :
*[*`t`,`y`] =**ode23s***(*`fun`,`trange`,`init`) - :
*[*`t`,`y`] =**ode23s***(*`fun`,`trange`,`init`,`ode_opt`) - :
*[*`t`,`y`] =**ode23s***(…,*`par1`,`par2`, …) - :
*[*`t`,`y`,`te`,`ye`,`ie`] =**ode23s***(…)* - :
`solution`=**ode23s***(…)* -
Solve a set of stiff Ordinary Differential Equations (stiff ODEs) with a Rosenbrock method of order (2,3).

`fun`is a function handle, inline function, or string containing the name of the function that defines the ODE:`M y' = f(t,y)`

. The function must accept two inputs where the first is time`t`and the second is a column vector of unknowns`y`.`M`is a constant mass matrix, non-singular and possibly sparse. Set the field`"Mass"`

in`odeopts`using`odeset`to specify a mass matrix.`trange`specifies the time interval over which the ODE will be evaluated. Typically, it is a two-element vector specifying the initial and final times (`[tinit, tfinal]`

). If there are more than two elements then the solution will also be evaluated at these intermediate time instances using an interpolation procedure of the same order as the one of the solver.By default,

`ode23s`

uses an adaptive timestep with the`integrate_adaptive`

algorithm. The tolerance for the timestep computation may be changed by using the options`"RelTol"`

and`"AbsTol"`

.`init`contains the initial value for the unknowns. If it is a row vector then the solution`y`will be a matrix in which each column is the solution for the corresponding initial value in`init`.The optional fourth argument

`ode_opt`specifies non-default options to the ODE solver. It is a structure generated by`odeset`

.`ode23s`

will ignore the following options:`"BDF"`

,`"InitialSlope"`

,`"MassSingular"`

,`"MStateDependence"`

,`"MvPattern"`

,`"MaxOrder"`

,`"Non-negative"`

.The function typically returns two outputs. Variable

`t`is a column vector and contains the times where the solution was found. The output`y`is a matrix in which each column refers to a different unknown of the problem and each row corresponds to a time in`t`. If`trange`specifies intermediate time steps, only those will be returned.The output can also be returned as a structure

`solution`which has a field`x`containing a row vector of times where the solution was evaluated and a field`y`containing the solution matrix such that each column corresponds to a time in`x`. Use`fieldnames (`

to see the other fields and additional information returned.`solution`)If using the

`"Events"`

option then three additional outputs may be returned.`te`holds the time when an Event function returned a zero.`ye`holds the value of the solution at time`te`.`ie`contains an index indicating which Event function was triggered in the case of multiple Event functions.Example: Solve the stiff Van der Pol equation

f = @(

`t`,`y`) [`y`(2); 1000*(1 -`y`(1)^2) *`y`(2) -`y`(1)]; opt = odeset ('Mass', [1 0; 0 1], 'MaxStep', 1e-1); [vt, vy] = ode23s (f, [0 2000], [2 0], opt);

- :
*[*`t`,`y`] =**ode15s***(*`fun`,`trange`,`y0`) - :
*[*`t`,`y`] =**ode15s***(*`fun`,`trange`,`y0`,`ode_opt`) - :
*[*`t`,`y`,`te`,`ye`,`ie`] =**ode15s***(…)* - :
`solution`=**ode15s***(…)* - :
**ode15s***(…)* Solve a set of stiff Ordinary Differential Equations (ODEs) or stiff semi-explicit index 1 Differential Algebraic Equations (DAEs).

`ode15s`

uses a variable step, variable order BDF (Backward Differentiation Formula) method that ranges from order 1 to 5.`fun`is a function handle, inline function, or string containing the name of the function that defines the ODE:`y' = f(t,y)`

. The function must accept two inputs where the first is time`t`and the second is a column vector of unknowns`y`.`trange`specifies the time interval over which the ODE will be evaluated. Typically, it is a two-element vector specifying the initial and final times (`[tinit, tfinal]`

). If there are more than two elements then the solution will also be evaluated at these intermediate time instances.`init`contains the initial value for the unknowns. If it is a row vector then the solution`y`will be a matrix in which each column is the solution for the corresponding initial value in`init`.The optional fourth argument

`ode_opt`specifies non-default options to the ODE solver. It is a structure generated by`odeset`

.The function typically returns two outputs. Variable

`t`is a column vector and contains the times where the solution was found. The output`y`is a matrix in which each column refers to a different unknown of the problem and each row corresponds to a time in`t`.`solution`which has a field`x`containing a row vector of times where the solution was evaluated and a field`y`containing the solution matrix such that each column corresponds to a time in`x`. Use`fieldnames (`

to see the other fields and additional information returned.`solution`)If no output arguments are requested, and no

`"OutputFcn"`

is specified in`ode_opt`, then the`"OutputFcn"`

is set to`odeplot`

and the results of the solver are plotted immediately.`"Events"`

option then three additional outputs may be returned.`te`holds the time when an Event function returned a zero.`ye`holds the value of the solution at time`te`.`ie`contains an index indicating which Event function was triggered in the case of multiple Event functions.Example: Solve Robertson’s equations:

function r = robertson_dae (

`t`,`y`) r = [ -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 ]; endfunction opt = odeset ("Mass", [1 0 0; 0 1 0; 0 0 0], "MStateDependence", "none"); [`t`,`y`] = ode15s (@robertson_dae, [0, 1e3], [1; 0; 0], opt);

- :
*[*`t`,`y`] =**ode15i***(*`fun`,`trange`,`y0`,`yp0`) - :
*[*`t`,`y`] =**ode15i***(*`fun`,`trange`,`y0`,`yp0`,`ode_opt`) - :
*[*`t`,`y`,`te`,`ye`,`ie`] =**ode15i***(…)* - :
`solution`=**ode15i***(…)* - :
**ode15i***(…)* Solve a set of fully-implicit Ordinary Differential Equations (ODEs) or index 1 Differential Algebraic Equations (DAEs).

`ode15i`

uses a variable step, variable order BDF (Backward Differentiation Formula) method that ranges from order 1 to 5.`fun`is a function handle, inline function, or string containing the name of the function that defines the ODE:`0 = f(t,y,yp)`

. The function must accept three inputs where the first is time`t`, the second is the function value`y`(a column vector), and the third is the derivative value`yp`(a column vector).`trange`specifies the time interval over which the ODE will be evaluated. Typically, it is a two-element vector specifying the initial and final times (`[tinit, tfinal]`

). If there are more than two elements then the solution will also be evaluated at these intermediate time instances.`y0`and`yp0`contain the initial values for the unknowns`y`and`yp`. If they are row vectors then the solution`y`will be a matrix in which each column is the solution for the corresponding initial value in`y0`and`yp0`.`y0`and`yp0`must be consistent initial conditions, meaning that`f(t,y0,yp0) = 0`

is satisfied. The function`decic`

may be used to compute consistent initial conditions given initial guesses.The optional fifth argument

`ode_opt`specifies non-default options to the ODE solver. It is a structure generated by`odeset`

.`t`is a column vector and contains the times where the solution was found. The output`y`is a matrix in which each column refers to a different unknown of the problem and each row corresponds to a time in`t`.`solution`which has a field`x`containing a row vector of times where the solution was evaluated and a field`y`containing the solution matrix such that each column corresponds to a time in`x`. Use`fieldnames (`

to see the other fields and additional information returned.`solution`)`"OutputFcn"`

is specified in`ode_opt`, then the`"OutputFcn"`

is set to`odeplot`

and the results of the solver are plotted immediately.`"Events"`

option then three additional outputs may be returned.`te`holds the time when an Event function returned a zero.`ye`holds the value of the solution at time`te`.`ie`contains an index indicating which Event function was triggered in the case of multiple Event functions.Example: Solve Robertson’s equations:

function r = robertson_dae (

`t`,`y`,`yp`) r = [ -(`yp`(1) + 0.04*`y`(1) - 1e4*`y`(2)*`y`(3)) -(`yp`(2) - 0.04*`y`(1) + 1e4*`y`(2)*`y`(3) + 3e7*`y`(2)^2)`y`(1) +`y`(2) +`y`(3) - 1 ]; endfunction [`t`,`y`] = ode15i (@robertson_dae, [0, 1e3], [1; 0; 0], [-1e-4; 1e-4; 0]);

- :
*[*`y0_new`,`yp0_new`] =**decic***(*`fun`,`t0`,`y0`,`fixed_y0`,`yp0`,`fixed_yp0`) - :
*[*`y0_new`,`yp0_new`] =**decic***(*`fun`,`t0`,`y0`,`fixed_y0`,`yp0`,`fixed_yp0`,`options`) - :
*[*`y0_new`,`yp0_new`,`resnorm`] =**decic***(…)* -
Compute consistent implicit ODE initial conditions

`y0_new`and`yp0_new`given initial guesses`y0`and`yp0`.A maximum of

`length (`

components between`y0`)`fixed_y0`and`fixed_yp0`may be chosen as fixed values.`fun`is a function handle. The function must accept three inputs where the first is time`t`, the second is a column vector of unknowns`y`, and the third is a column vector of unknowns`yp`.`t0`is the initial time such that

, specified as a scalar.`fun`(`t0`,`y0_new`,`yp0_new`) = 0`y0`is a vector used as the initial guess for`y`.`fixed_y0`is a vector which specifies the components of`y0`to hold fixed. Choose a maximum of`length (`

components between`y0`)`fixed_y0`and`fixed_yp0`as fixed values. Set`fixed_y0`(i) component to 1 if you want to fix the value of`y0`(i). Set`fixed_y0`(i) component to 0 if you want to allow the value of`y0`(i) to change.`yp0`is a vector used as the initial guess for`yp`.`fixed_yp0`is a vector which specifies the components of`yp0`to hold fixed. Choose a maximum of`length (`

components between`yp0`)`fixed_y0`and`fixed_yp0`as fixed values. Set`fixed_yp0`(i) component to 1 if you want to fix the value of`yp0`(i). Set`fixed_yp0`(i) component to 0 if you want to allow the value of`yp0`(i) to change.The optional seventh argument

`options`is a structure array. Use`odeset`

to generate this structure. The relevant options are`RelTol`

and`AbsTol`

which specify the error thresholds used to compute the initial conditions.The function typically returns two outputs. Variable

`y0_new`is a column vector and contains the consistent initial value of`y`. The output`yp0_new`is a column vector and contains the consistent initial value of`yp`.The optional third output

`resnorm`is the norm of the vector of residuals. If`resnorm`is small,`decic`

has successfully computed the initial conditions. If the value of`resnorm`is large, use`RelTol`

and`AbsTol`

to adjust it.Example: Compute initial conditions for Robertson’s equations:

function r = robertson_dae (

`t`,`y`,`yp`) r = [ -(`yp`(1) + 0.04*`y`(1) - 1e4*`y`(2)*`y`(3)) -(`yp`(2) - 0.04*`y`(1) + 1e4*`y`(2)*`y`(3) + 3e7*`y`(2)^2)`y`(1) +`y`(2) +`y`(3) - 1 ]; endfunction[

`y0_new`,`yp0_new`] = decic (@robertson_dae, 0, [1; 0; 0], [1; 1; 0], [-1e-4; 1; 0], [0; 0; 0]);

- :
`odestruct`=**odeset***()* - :
`odestruct`=**odeset***(*`"field1"`,`value1`,`"field2"`,`value2`, …) - :
`odestruct`=**odeset***(*`oldstruct`,`"field1"`,`value1`,`"field2"`,`value2`, …) - :
`odestruct`=**odeset***(*`oldstruct`,`newstruct`) - :
**odeset***()* -
Create or modify an ODE options structure.

When called with no input argument and one output argument, return a new ODE options structure that contains all possible fields initialized to their default values. If no output argument is requested, display a list of the common ODE solver options along with their default value.

If called with name-value input argument pairs

`"field1"`,`"value1"`,`"field2"`,`"value2"`, … return a new ODE options structure with all the most common option fields initialized,**and**set the values of the fields`"field1"`,`"field2"`, … to the values`value1`,`value2`, ….If called with an input structure

`oldstruct`then overwrite the values of the options`"field1"`,`"field2"`, … with new values`value1`,`value2`, … and return the modified structure.When called with two input ODE options structures

`oldstruct`and`newstruct`overwrite all values from the structure`oldstruct`with new values from the structure`newstruct`. Empty values in`newstruct`will not overwrite values in`oldstruct`.The most commonly used ODE options, which are always assigned a value by

`odeset`

, are the following:`AbsTol`

: positive scalar | vector, def.`1e-6`

Absolute error tolerance.

`BDF`

: {`"off"`

} |`"on"`

Use BDF formulas in implicit multistep methods.

*Note*: This option is not yet implemented.`Events`

: function_handleEvent function. An event function must have the form

`[value, isterminal, direction] = my_events_f (t, y)`

`InitialSlope`

: vectorConsistent initial slope vector for DAE solvers.

`InitialStep`

: positive scalarInitial time step size.

`Jacobian`

: matrix | function_handleJacobian matrix, specified as a constant matrix or a function of time and state.

`JConstant`

: {`"off"`

} |`"on"`

Specify whether the Jacobian is a constant matrix or depends on the state.

`JPattern`

: sparse matrixIf the Jacobian matrix is sparse and non-constant but maintains a constant sparsity pattern, specify the sparsity pattern.

`Mass`

: matrix | function_handleMass matrix, specified as a constant matrix or a function of time and state.

`MassSingular`

: {`"maybe"`

} |`"yes"`

|`"on"`

Specify whether the mass matrix is singular.

`MaxOrder`

: {`5`

} |`4`

|`3`

|`2`

|`1`

Maximum order of formula.

`MaxStep`

: positive scalarMaximum time step value.

`MStateDependence`

: {`"weak"`

} |`"none"`

|`"strong"`

Specify whether the mass matrix depends on the state or only on time.

`MvPattern`

: sparse matrixIf the mass matrix is sparse and non-constant but maintains a constant sparsity pattern, specify the sparsity pattern.

*Note*: This option is not yet implemented.`NonNegative`

: scalar | vectorSpecify elements of the state vector that are expected to remain non-negative during the simulation.

`NormControl`

: {`"off"`

} |`"on"`

Control error relative to the 2-norm of the solution, rather than its absolute value.

`OutputFcn`

: function_handleFunction to monitor the state during the simulation. For the form of the function to use see

`odeplot`

.`OutputSel`

: scalar | vectorIndices of elements of the state vector to be passed to the output monitoring function.

`Refine`

: positive scalarSpecify whether output should be returned only at the end of each time step or also at intermediate time instances. The value should be a scalar indicating the number of equally spaced time points to use within each timestep at which to return output.

*Note*: This option is not yet implemented.`RelTol`

: positive scalarRelative error tolerance.

`Stats`

: {`"off"`

} |`"on"`

Print solver statistics after simulation.

`Vectorized`

: {`"off"`

} |`"on"`

Specify whether

`odefun`

can be passed multiple values of the state at once.

Field names that are not in the above list are also accepted and added to the result structure.

**See also:**odeget.

- :
`val`=**odeget***(*`ode_opt`,`field`) - :
`val`=**odeget***(*`ode_opt`,`field`,`default`) -
Query the value of the property

`field`in the ODE options structure`ode_opt`.If called with two input arguments and the first input argument

`ode_opt`is an ODE option structure and the second input argument`field`is a string specifying an option name, then return the option value`val`corresponding to`field`from`ode_opt`.If called with an optional third input argument, and

`field`is not set in the structure`ode_opt`, then return the default value`default`instead.**See also:**odeset.

- :
`stop_solve`=**odeplot***(*`t`,`y`,`flag`) -
Open a new figure window and plot the solution of an ode problem at each time step during the integration.

The types and values of the input parameters

`t`and`y`depend on the input`flag`that is of type string. Valid values of`flag`are:`"init"`

The input

`t`must be a column vector of length 2 with the first and last time step (`[`

. The input`tfirst``tlast`]`y`contains the initial conditions for the ode problem (`y0`).`""`

The input

`t`must be a scalar double specifying the time for which the solution in input`y`was calculated.`"done"`

The inputs should be empty, but are ignored if they are present.

`odeplot`

always returns false, i.e., don’t stop the ode solver.Example: solve an anonymous implementation of the

`"Van der Pol"`

equation and display the results while solving.fvdp = @(t,y) [y(2); (1 - y(1)^2) * y(2) - y(1)]; opt = odeset ("OutputFcn", @odeplot, "RelTol", 1e-6); sol = ode45 (fvdp, [0 20], [2 0], opt);

Background Information: This function is called by an ode solver function if it was specified in the

`"OutputFcn"`

property of an options structure created with`odeset`

. The ode solver will initially call the function with the syntax`odeplot ([`

. The function initializes internal variables, creates a new figure window, and sets the x limits of the plot. Subsequently, at each time step during the integration the ode solver calls`tfirst`,`tlast`],`y0`, "init")`odeplot (`

. At the end of the solution the ode solver calls`t`,`y`, [])`odeplot ([], [], "done")`

so that odeplot can perform any clean-up actions required.

Previous: Differential-Algebraic Equations, Up: Differential Equations [Contents][Index]