Matlab-compatible solvers (GNU Octave (version 9.1.0)) (2024)

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

24.3 Matlab-compatible solvers ¶

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

odeset
odeget

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.
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 andM. W. Reichelt, The MATLAB ODE Suite, SIAM Journal onScientific Computing, Vol. 18, 1997, pp. 1–22,DOI: https://doi.org/10.1137/S1064827594276424.

: [t, y] = ode45 (fcn, trange, init) ¶

: [t, y] = ode45 (fcn, 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.

fcn is a function handle, inline function, or string containing thename of the function that defines the ODE: y' = f(t,y). The functionmust accept two inputs where the first is time t and the second is acolumn vector of unknowns y.

trange specifies the time interval over which the ODE will beevaluated. Typically, it is a two-element vector specifying the initial andfinal times ([tinit, tfinal]). If there are more than two elementsthen the solution will also be evaluated at these intermediate timeinstances.

By default, ode45 uses an adaptive timestep with theintegrate_adaptive algorithm. The tolerance for the timestepcomputation may be changed by using the options "RelTol" and"AbsTol".

init contains the initial value for the unknowns. If it is a rowvector then the solution y will be a matrix in which each column isthe solution for the corresponding initial value in init.

The optional fourth argument ode_opt specifies non-default options tothe ODE solver. It is a structure generated by odeset.

The function typically returns two outputs. Variable t is acolumn vector and contains the times where the solution was found. Theoutput y is a matrix in which each column refers to a differentunknown of the problem and each row corresponds to a time in t.

The output can also be returned as a structure solution which has afield x containing a row vector of times where the solution wasevaluated and a field y containing the solution matrix such that eachcolumn corresponds to a time in x. Usefieldnames(solution) to see the other fields andadditional information returned.

If no output arguments are requested, and no "OutputFcn" isspecified in ode_opt, then the "OutputFcn" is set toodeplot and the results of the solver are plotted immediately.

If using the "Events" option then three additional outputs may bereturned. te holds the time when an Event function returned a zero.ye holds the value of the solution at time te. iecontains an index indicating which Event function was triggered in the caseof 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]);

See also: odeset, odeget, ode23, ode15s.

: [t, y] = ode23 (fcn, trange, init) ¶

: [t, y] = ode23 (fcn, 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.

By default, ode23 uses an adaptive timestep with theintegrate_adaptive algorithm. The tolerance for the timestepcomputation may be changed by using the options "RelTol" and"AbsTol".

init contains the initial value for the unknowns. If it is a rowvector then the solution y will be a matrix in which each column isthe solution for the corresponding initial value in init.

The optional fourth argument ode_opt specifies non-default options tothe ODE solver. It is a structure generated by odeset.

If no output arguments are requested, and no "OutputFcn" isspecified in ode_opt, then the "OutputFcn" is set toodeplot and the results of the solver are plotted immediately.

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 seehttps://en.wikipedia.org/wiki/List_of_Runge%E2%80%93Kutta_methods.

See also: odeset, odeget, ode45, ode15s.

: [t, y] = ode23s (fcn, trange, init) ¶

: [t, y] = ode23s (fcn, 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 aRosenbrock method of order (2,3).

fcn is a function handle, inline function, or string containing thename of the function that defines the ODE: M y' = f(t,y). Thefunction must accept two inputs where the first is time t and thesecond is a column vector of unknowns y. M is a constant massmatrix, non-singular and possibly sparse. Set the field "Mass" inodeopts using odeset to specify a mass matrix.

trange specifies the time interval over which the ODE will beevaluated. Typically, it is a two-element vector specifying the initialand final times ([tinit, tfinal]). If there are more than twoelements then the solution will also be evaluated at these intermediatetime instances using an interpolation procedure of the same order as theone of the solver.

By default, ode23s uses an adaptive timestep with theintegrate_adaptive algorithm. The tolerance for the timestepcomputation may be changed by using the options "RelTol" and"AbsTol".

init contains the initial value for the unknowns. If it is a rowvector then the solution y will be a matrix in which each column isthe solution for the corresponding initial value in init.

The optional fourth argument ode_opt specifies non-default options tothe ODE solver. It is a structure generated by odeset.ode23s will ignore the following options: "BDF","InitialSlope", "MassSingular", "MStateDependence","MvPattern", "MaxOrder", "Non-negative".

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);

See also: odeset, daspk, dassl.

: [t, y] = ode15s (fcn, trange, y0) ¶

: [t, y] = ode15s (fcn, trange, y0, ode_opt) ¶

: [t, y, te, ye, ie] = ode15s (…) ¶

: solution = ode15s (…) ¶

: ode15s (…) ¶

Solve a set of stiff Ordinary Differential Equations (ODEs) or stiffsemi-explicit index 1 Differential Algebraic Equations (DAEs).

ode15s uses a variable step, variable order BDF (BackwardDifferentiation Formula) method that ranges from order 1 to 5.

init contains the initial value for the unknowns. If it is a rowvector then the solution y will be a matrix in which each column isthe solution for the corresponding initial value in init.

The optional fourth argument ode_opt specifies non-default options tothe ODE solver. It is a structure generated by odeset.

If no output arguments are requested, and no "OutputFcn" isspecified in ode_opt, then the "OutputFcn" is set toodeplot and the results of the solver are plotted immediately.

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)^2y(1) + y(2) + y(3) - 1 ];endfunctionopt = odeset ("Mass", [1 0 0; 0 1 0; 0 0 0], "MStateDependence", "none");[t,y] = ode15s (@robertson_dae, [0, 1e3], [1; 0; 0], opt);

See also: decic, odeset, odeget, ode23, ode45.

: [t, y] = ode15i (fcn, trange, y0, yp0) ¶

: [t, y] = ode15i (fcn, trange, y0, yp0, ode_opt) ¶

: [t, y, te, ye, ie] = ode15i (…) ¶

: solution = ode15i (…) ¶

: ode15i (…) ¶

Solve a set of fully-implicit Ordinary Differential Equations (ODEs) orindex 1 Differential Algebraic Equations (DAEs).

ode15i uses a variable step, variable order BDF (BackwardDifferentiation Formula) method that ranges from order 1 to 5.

fcn is a function handle, inline function, or string containing thename of the function that defines the ODE: 0 = f(t,y,yp). Thefunction must accept three inputs where the first is time t, thesecond is the function value y (a column vector), and the thirdis the derivative value yp (a column vector).

y0 and yp0 contain the initial values for the unknowns yand yp. If they are row vectors then the solution y will be amatrix in which each column is the solution for the corresponding initialvalue in y0 and yp0.

y0 and yp0 must be consistent initial conditions, meaning thatf(t,y0,yp0) = 0 is satisfied. The function decic may be usedto compute consistent initial conditions given initial guesses.

The optional fifth argument ode_opt specifies non-default options tothe ODE solver. It is a structure generated by odeset.

If no output arguments are requested, and no "OutputFcn" isspecified in ode_opt, then the "OutputFcn" is set toodeplot and the results of the solver are plotted immediately.

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]);

See also: decic, odeset, odeget.

: [y0_new, yp0_new] = decic (fcn, t0, y0, fixed_y0, yp0, fixed_yp0) ¶

: [y0_new, yp0_new] = decic (fcn, t0, y0, fixed_y0, yp0, fixed_yp0, options) ¶

: [y0_new, yp0_new, resnorm] = decic (…) ¶

Compute consistent implicit ODE initial conditions y0_new andyp0_new given initial guesses y0 and yp0.

A maximum of length (y0) components between fixed_y0 andfixed_yp0 may be chosen as fixed values.

fcn is a function handle. The function must accept three inputs wherethe first is time t, the second is a column vector of unknownsy, and the third is a column vector of unknowns yp.

t0 is the initial time such thatfcn(t0, y0_new, yp0_new) = 0, specified as ascalar.

y0 is a vector used as the initial guess for y.

fixed_y0 is a vector which specifies the components of y0 tohold fixed. Choose a maximum of length (y0) components betweenfixed_y0 and fixed_yp0 as fixed values.Set fixed_y0(i) component to 1 if you want to fix the value ofy0(i).Set fixed_y0(i) component to 0 if you want to allow the value ofy0(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 tohold fixed. Choose a maximum of length (yp0) componentsbetween fixed_y0 and fixed_yp0 as fixed values.Set fixed_yp0(i) component to 1 if you want to fix the value ofyp0(i).Set fixed_yp0(i) component to 0 if you want to allow the value ofyp0(i) to change.

The optional seventh argument options is a structure array. Useodeset to generate this structure. The relevant options areRelTol and AbsTol which specify the error thresholds used tocompute the initial conditions.

The function typically returns two outputs. Variable y0_new is acolumn vector and contains the consistent initial value of y. Theoutput yp0_new is a column vector and contains the consistent initialvalue of yp.

The optional third output resnorm is the norm of the vector ofresiduals. If resnorm is small, decic has successfullycomputed 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]);

See also: ode15i, odeset.

: 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 ODEoptions structure that contains all possible fields initialized to theirdefault values. If no output argument is requested, display a list ofthe common ODE solver options along with their default value.

If called with name-value input argument pairs "field1","value1", "field2", "value2", … return a newODE options structure with all the most common option fieldsinitialized, and set the values of the fields "field1","field2", … to the values value1, value2,….

If called with an input structure oldstruct then overwrite thevalues of the options "field1", "field2", … withnew values value1, value2, … and return themodified structure.

When called with two input ODE options structures oldstruct andnewstruct overwrite all values from the structureoldstruct with new values from the structure newstruct.Empty values in newstruct will not overwrite values inoldstruct.

The most commonly used ODE options, which are always assigned a valueby 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_handle: Event function. An event function must have the form[value, isterminal, direction] = my_events_f (t, y)
InitialSlope: vector: Consistent initial slope vector for DAE solvers.
InitialStep: positive scalar: Initial time step size.
Jacobian: matrix | function_handle: Jacobian matrix, specified as a constant matrix or a function oftime and state.
JConstant: {"off"} | "on": Specify whether the Jacobian is a constant matrix or depends on thestate.
JPattern: sparse matrix: If the Jacobian matrix is sparse and non-constant but maintains aconstant sparsity pattern, specify the sparsity pattern.
Mass: matrix | function_handle: Mass matrix, specified as a constant matrix or a function oftime 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 scalar: Maximum time step value.
MStateDependence: {"weak"} | "none" | "strong": Specify whether the mass matrix depends on the state or only on time.
MvPattern: sparse matrix: If the mass matrix is sparse and non-constant but maintains aconstant sparsity pattern, specify the sparsity pattern.Note: This option is not yet implemented.
NonNegative: scalar | vector: Specify elements of the state vector that are expected to remainnon-negative during the simulation.
NormControl: {"off"} | "on": Control error relative to the 2-norm of the solution, rather than itsabsolute value.
OutputFcn: function_handle: Function to monitor the state during the simulation. For the form ofthe function to use see odeplot.
OutputSel: scalar | vector: Indices of elements of the state vector to be passed to the outputmonitoring function.
Refine: positive scalar: Specify whether output should be returned only at the end of eachtime step or also at intermediate time instances. The value should bea scalar indicating the number of equally spaced time points to usewithin each timestep at which to return output.
RelTol: positive scalar: Relative error tolerance.
Stats: {"off"} | "on": Print solver statistics after simulation.
Vectorized: {"off"} | "on": Specify whether odefcn can be passed multiple values of thestate at once.

Field names that are not in the above list are also accepted andadded 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 structureode_opt.

If called with two input arguments and the first input argumentode_opt is an ODE option structure and the second input argumentfield is a string specifying an option name, then return the optionvalue val corresponding to field from ode_opt.

If called with an optional third input argument, and field isnot set in the structure ode_opt, then return the default valuedefault instead.

See also: odeset.

: stop_solve = odeplot (t, y, flag) ¶

Open a new figure window and plot the solution of an ode problem at eachtime step during the integration.

The types and values of the input parameters t and y depend onthe input flag that is of type string. Valid values of flagare:

"init": The input t must be a column vector of length 2 with the first andlast time step ([tfirst tlast]. The input ycontains the initial conditions for the ode problem (y0).
"": The input t must be a scalar double or vector specifying the time(s)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 whilesolving.

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 inthe "OutputFcn" property of an options structure created withodeset. The ode solver will initially call the function with thesyntax odeplot ([tfirst, tlast], y0, "init"). Thefunction initializes internal variables, creates a new figure window, andsets the x limits of the plot. Subsequently, at each time step during theintegration the ode solver calls odeplot (t, y, []).At the end of the solution the ode solver callsodeplot ([], [], "done") so that odeplot can perform any clean-upactions required.

See also: odeset, odeget, ode23, ode45.