# ode

Ordinary differential equations

Since R2023b

## Description

An `ode` object defines a system of ordinary differential equations or differential algebraic equations to solve.

You can solve initial value problems of the form $y\text{'}=f\left(t,y\right)$, $f\left(t,y,y\text{'}\right)=0$, or problems that involve a mass matrix, $M\left(t,y\right)y\text{'}=f\left(t,y\right)$.

Define aspects of the problem using properties of the `ode` object, such as `ODEFcn`, `InitialTime`, and `InitialValue`. You can select a specific solver to use, or let MATLAB® choose an appropriate solver based on properties of the equations. After you create an `ode` object, you can solve the equations using the `solve` or `solutionFcn` object functions.

## Creation

### Syntax

``F = ode``
``F = ode(Name=Value)``

### Description

````F = ode` creates an `ode` object with default properties.```

example

````F = ode(Name=Value)` specifies one or more property values using name-value arguments. For example, you can specify the equations to be solved, the initial time for integration, and the value of the solution at the initial time using the `ODEFcn`, `InitialTime`, and `InitialValue` properties.```

example

## Properties

expand all

### Problem Definition

Equations to solve, specified as a function handle that defines the system of differential equations to be integrated. The function handle can be an anonymous function or a handle to a named function file.

• The function `dydt = ODEFcn(t,y)`, for a scalar `t` and a column vector `y`, must return a column vector `dydt` that corresponds to $f\left(t,y\right)$. `ODEFcn` must accept at least two input arguments, `t` and `y`, even if one of the arguments is not used in the function.

• The function ```f = ODEFcn(t,y,yp)```, when `EquationType` is `"fullyimplicit"` and for a scalar `t` and column vectors `y` and `yp`, must return a column vector `f` that corresponds to $f\left(t,y,y\text{'}\right)$. `ODEFcn` must accept at least three input arguments, `t`, `y`, and `yp`, even if one of the arguments is not used in the function. (since R2024b)

To supply parameter values, define `ODEFcn` as a function that accepts three inputs, `dydt = ODEFcn(t,y,p)`, or four inputs, `f = ODEFcn(t,y,yp,p)`, and then use the `Parameters` property to store parameter values (such as in a `struct` or `cell` array).

For example, to solve a single equation $y\text{'}=5y-3$, you can use the anonymous function ```@(t,y) 5*y-3```.

For a system of equations, the output of `ODEFcn` is a vector. Each element in the vector is the computed value of the right side of one equation. For example, consider this system of two equations.

`$\begin{array}{l}y{\text{'}}_{1}={y}_{1}+2{y}_{2}\\ y{\text{'}}_{2}=3{y}_{1}+2{y}_{2}\end{array}$`

An anonymous function that defines these equations is ```@(t,y) [y(1)+2*y(2); 3*y(1)+2*y(2)]```.

Example: `F.ODEFcn = @myFcn` specifies a handle to a named function file `myFcn.m` containing the equations.

Example: `F.ODEFcn = @(t,y) [y(2); -y(1)]` specifies an anonymous function that defines a system of two equations.

Example: `F.ODEFcn = @(t,y,yp) [yp(1)-y(2); yp(2)+1]` specifies an anonymous function that defines a system of two equations.

Example: `F.ODEFcn = @(t,y,p) 5*y*p(1)-3*p(2)` specifies an anonymous function for a single equation that uses two parameters.

Data Types: `function_handle`

Initial time for integration, specified as a real scalar. The value of `InitialTime` is the beginning of the integration interval where the initial conditions specified in `InitialValue` are applied by the solver before beginning integration steps.

Example: `F.InitialTime = 10;`

Data Types: `single` | `double`

Value of solution at `InitialTime`, specified as a scalar or vector. `InitialValue` must be a vector with the same length as the output of `ODEFcn` so that an initial value is specified for each equation defined in `ODEFcn`. The value of `InitialTime` is the beginning of the integration interval where the initial conditions specified in `InitialValue` are applied by the solver before beginning integration steps.

Example: `F.InitialValue = 1;`

Example: `F.InitialValue = [1 2 3];`

Data Types: `single` | `double`

Equation parameters, specified as an array of any size or data type. The values you store in `Parameters` can be supplied to any of the functions used for the `ODEFcn`, `Jacobian`, `MassMatrix`, or `EventDefinition` properties by specifying an extra input argument in the function.

For instance, you can supply parameter values stored in the `Parameters` property to `ODEFcn` by specifying `ODEFcn` as a function handle that accepts three inputs, `dydt = ODEFcn(t,y,p)`. For example, ```F.ODEFcn = @(t,y,p) 5*y*p(1)-3*p(2)``` specifies an anonymous function for a single equation that uses two parameters. If you then specify `F.Parameters = [2 3]`, then `ODEFcn` uses the parameter values `p(1) = 2` and `p(2) = 3` whenever the solver calls the function.

Example: `F.Parameters = [0.1 0.5]` specifies two parameter values.

Mass matrix, specified as an `odeMassMatrix` object, matrix, or handle to a function that evaluates the mass matrix.

`ode` objects can represent problems of the form $M\left(t,y\right)\text{\hspace{0.17em}}y\text{'}=f\left(t,y\right)$, where $M\left(t,y\right)$ is a mass matrix that can be full or sparse. The mass matrix encodes linear combinations of derivatives on the left side of the equation.

• When the mass matrix is nonsingular, the equation simplifies to $y\text{'}={M}^{-1}\text{\hspace{0.17em}}f\left(t,y\right)$ and the ODE has a solution for any initial value. However, it is often more convenient and natural to express the model in terms of the mass matrix directly using $M\left(t,y\right)\text{\hspace{0.17em}}y\text{'}=f\left(t,y\right)$, and avoiding the computation of the matrix inverse reduces the storage and execution time needed to solve the problem.

• When the mass matrix is singular, then the problem is a system of differential algebraic equations (DAEs). A DAE has a solution only when the initial values are consistent; that is, you must specify the initial slope $y{\text{'}}_{0}$ using the `InitialSlope` property such that the initial conditions are all consistent, $M\left({t}_{0},{y}_{0}\right)y{\text{'}}_{0}=f\left({t}_{0},{y}_{0}\right)$. If the specified initial conditions are not consistent with the `InitialTime`, `InitialValue`, and `MassMatrix` properties, then the solver treats them as guesses and attempts to compute consistent values for the initial slopes that are close to the guesses before continuing to solve the problem. For more information, see Solve Differential Algebraic Equations (DAEs).

You can specify the value of the `MassMatrix` property as:

• An `odeMassMatrix` object, which represents the mass matrix and its associated properties. You can specify whether the mass matrix is singular or has state dependence.

• A constant matrix with calculated values for $M\left(t,y\right)$.

• A handle to a function that computes the matrix elements and that accepts two input arguments, `M = Mass(t,y)`. To give the function access to parameter values stored in the `Parameters` property, specify a third input argument in the function definition, ```M = Mass(t,y,p)```.

If you specify a matrix or function handle, then MATLAB converts it to an `odeMassMatrix` object.

Example: `F.MassMatrix = @Mass` specifies the function `Mass` that evaluates the mass matrix.

Example: `F.MassMatrix = [1 0; -2 1]` specifies a constant mass matrix.

Example: ```F.MassMatrix = odeMassMatrix(MassMatrix=@Mass,StateDependence="strong",SparsityPattern=S)``` specifies a state-dependent mass matrix and the sparsity pattern when the mass matrix multiplies a vector.

Since R2024b

Equation form, specified as `"standard"` by default or `"fullyimplicit"` for implicit ODEs, that is, a system of differential equations of the form $f\left(t,y,y\text{'}\right)=0$. If you specify the Jacobian when the `EquationType` property is `"fullyimplicit"`, the Jacobian must be a function handle or cell array.

Jacobian matrix, specified as an `odeJacobian` object, matrix, cell array, or handle to a function that evaluates the Jacobian. The Jacobian is a matrix of partial derivatives of the functions that define the system of differential equations.

`$J=\frac{\partial f}{\partial y}=\left[\begin{array}{ccc}\frac{\partial {f}_{1}}{\partial {y}_{1}}& \frac{\partial {f}_{1}}{\partial {y}_{2}}& \cdots \\ \frac{\partial {f}_{2}}{\partial {y}_{1}}& \frac{\partial {f}_{2}}{\partial {y}_{2}}& \cdots \\ ⋮& ⋮& \end{array}\right]$`

For stiff ODE problems, providing information about the Jacobian matrix is critical for reliability and efficiency of the solver. If you do not provide the Jacobian, then the ODE solver approximates it numerically using finite differences.

For large systems of equations where it is not feasible to provide the entire analytic Jacobian, you can specify the sparsity pattern of the Jacobian matrix instead. The solver uses the sparsity pattern to calculate a sparse Jacobian.

You can specify the `Jacobian` property as:

• An `odeJacobian` object, which can represent either the Jacobian matrix or its sparsity pattern.

• A constant matrix with calculated values for $\frac{\partial f}{\partial y}$.

• When `EquationType` is `"fullyimplicit"`, a two-element cell array with calculated values for the constant Jacobian with respect to `y` in the first element and `yp` in the second element. If you specify one of the elements as `[]`, the ODE solver approximates the corresponding Jacobian numerically while taking the provided values in the other element into account. (since R2024b)

• A handle to a function that computes the matrix elements and that accepts two input arguments, `dfdy = Fjac(t,y)`. To give the function access to parameter values in the `Parameters` property, specify a third input argument in the function definition, ```dfdy = Fjac(t,y,p)```.

• When `EquationType` is `"fullyimplicit"`, a handle to a function that computes the matrix elements and that accepts three input arguments, ```[dfdy,dfdp] = Fjac(t,y,yp)```. To give the function access to parameter values in the `Parameters` property, specify a fourth input argument in the function definition, `[dfdy,dfdp] = Fjac(t,y,yp,p)`. (since R2024b)

If you specify a matrix, function handle, or cell array, then MATLAB converts it to an `odeJacobian` object.

Example: `F.Jacobian = @Fjac` specifies the function `Fjac` that evaluates the Jacobian matrix.

Example: `F.Jacobian = [0 1; -2 1]` specifies a constant Jacobian matrix.

Example: `F.Jacobian = odeJacobian(SparsityPattern=S)` specifies the Jacobian sparsity pattern using sparse matrix `S`.

Events to detect, specified as an `odeEvent` object. Create an `odeEvent` object to define expressions that trigger an event when they cross zero. You can specify the direction of the zero crossing and what to do when an event triggers, including the use of a callback function.

Nonnegative solution components, specified as a scalar index or vector of indices. Use the `NonNegativeVariables` property to specify which solutions the solver must keep nonnegative. If `dydt = ODEFcn(t,y)`, then the indices correspond to elements in the vector `y`. For example, if you specify a value of `3`, then the solver keeps the solution component `y(3)` nonnegative.

Note

`NonNegativeVariables` is not available for the `ode23s` solver. Additionally, for `ode15s`, `ode23t`, and `ode23tb` the property is not available for problems that have a mass matrix.

Example: `F.NonNegativeVariables = [1 3]` specifies that the first and third solution components must be kept nonnegative.

Data Types: `single` | `double` | `int8` | `int16` | `int32` | `int64` | `uint8` | `uint16` | `uint32` | `uint64`

Initial value of `dy/dt`, specified as a vector. Use this property with the `ode15s` and `ode23t` solvers when solving DAEs. The specified vector is the initial slope $y{\text{'}}_{0}$ such that $M\left({t}_{0},{y}_{0}\right)y{\text{'}}_{0}=f\left({t}_{0},{y}_{0}\right)$. If the specified initial conditions are not consistent with the `InitialTime`, `InitialValue`, and `MassMatrix` properties, then the solver treats them as guesses and attempts to compute consistent values for the initial slopes that are close to the guesses before continuing to solve the problem. For more information, see Solve Differential Algebraic Equations (DAEs).

Data Types: `single` | `double`

Sensitivity analysis information, specified as an `odeSensitivity` object. The ODE solver performs sensitivity analysis only when the `Sensitivity` property is nonempty. The analysis is with respect to the `Parameters` property, which must be nonempty and numeric. See `odeSensitivity` for more information.

Example: `F.Sensitivity = odeSensitivity` performs sensitivity analysis on all parameters.

Example: `F.Sensitivity = odeSensitivity(ParameterIndices=1:2)` performs sensitivity analysis on the first two parameters.

### Solver Control

Solver method, specified as one of the values in this table. To change the solver based on specified parameters after creating an `ode` object, use the `selectSolver` function. For example, you can incorporate a stiffness detection heuristic when selecting a solver by using the `selectSolver` function with its `DetectStiffness` name-value argument.

ValueDescription

`"auto"`

Automatically selects a solver based on available problem information and tolerances. The `SelectedSolver` property shows the currently chosen solver, and the solver choice can change when you update property values of the `ode` object.

`"stiff"`

Similar to `"auto"` but chooses only from the stiff solvers: `"ode15s"`, `"ode23s"`, `"ode23t"`, `"ode23tb"`, `"cvodesstiff"`, `"idas"`.

`"nonstiff"`

Similar to `"auto"` but chooses only from the nonstiff solvers: `"ode23"`, `"ode45"`, `"ode78"`, `"ode89"`, `"ode113"`, `"cvodesnonstiff"`.

`"ode23"`

Runge-Kutta (2,3) pair. See `ode23`.

`"ode45"`

Runge-Kutta (4,5) pair. See `ode45`.

`"ode78"`

Runge-Kutta 8(7) pair with a 7th-order continuous extension. See `ode78`.

`"ode89"`

Runge-Kutta 9(8) pair with an 8th-order continuous extension. See `ode89`.

`"ode113"`

Variable-step, variable-order (VSVO) solver of orders 1 to 13. See `ode113`.

`"ode15s"`

Variable-step, variable-order (VSVO) solver based on the numerical differentiation formulas (NDFs) of orders 1 to 5. See `ode15s`.

`"ode23s"`

Modified Rosenbrock formula of order 2. See `ode23s`.

`"ode23t"`

Trapezoidal rule using a “free” interpolant. See `ode23t`.

`"ode23tb"`

Implicit Runge-Kutta formula with a trapezoidal rule step as its first stage and a backward differentiation formula of order 2 as its second stage. See `ode23tb`.

`"ode15i"` (since R2024b)

Variable-step, variable order (VSVO) solver based on the backward differentiation formulas (BDFs) of orders 1 to 5. See `ode15i`.

`"cvodesnonstiff"`

Variable-step, variable-order (VSVO) solver using Adams-Moulton formulas, with the order varying between 1 and 12. Supports sensitivity analysis. See CVODE and CVODES.

`"cvodesstiff"`

Variable-step, variable-order (VSVO) solver using backward differentiation formulas (BDFs) in fixed-leading coefficient form, with order varying between 1 and 5. Supports sensitivity analysis. See CVODE and CVODES.

`"idas"`

Variable-order, variable-coefficient solver using backward differentiation formulas (BDFs) in fixed-leading coefficient form, with order varying between 1 and 5. Supports sensitivity analysis. See IDA and IDAS.

Selected solver, returned as the name of a solver. If the value of `Solver` is `"auto"`, `"stiff"`, or `"nonstiff"`, then the `ode` object automatically chooses a solver and sets the value of `SelectedSolver` to the name of the chosen solver. If you manually select a solver, then `SelectedSolver` and `Solver` have the same value.

Solver-specific options, specified as a `matlab.ode.options.*` object. The `ode` object automatically populates the `SolverOptions` property with an options object appropriate for the selected solver. You can check available options for an existing `ode` object with the command `properties(F.SolverOptions)`. Specify options for the ODE problem by changing property values of the `matlab.ode.options.*` object. For example, `F.SolverOptions.OutputFcn = @odeplot` specifies an output function that the solver calls after each successful time step.

When the `Solver` property is `"auto"`, `"stiff"`, or `"nonstiff"`, the `SolverOptions` property is read-only.

This table summarizes the available options for each solver.

Options ObjectProperties

`matlab.ode.options.ODE23 Properties`

`matlab.ode.options.ODE45 Properties`

`matlab.ode.options.ODE78 Properties`

`matlab.ode.options.ODE89 Properties`

`matlab.ode.options.ODE113 Properties`

• `InitialStep`

• `MaxStep`

• `MinStep`

• `NormControl`

• `OutputFcn`

• `OutputSelection`

`matlab.ode.options.ODE15s Properties`

• `InitialStep`

• `MaxStep`

• `MinStep`

• `NormControl`

• `OutputFcn`

• `OutputSelection`

• `Vectorization`

• `BDF`

• `MaxOrder`

`matlab.ode.options.ODE15i Properties` (since R2024b)

• `InitialStep`

• `MaxStep`

• `MinStep`

• `NormControl`

• `OutputFcn`

• `OutputSelection`

• `Vectorization`

• `MaxOrder`

• `ComputeConsistentInitialConditions`

`matlab.ode.options.ODE23s Properties`

`matlab.ode.options.ODE23t Properties`

`matlab.ode.options.ODE23tb Properties`

• `InitialStep`

• `MaxStep`

• `MinStep`

• `NormControl`

• `OutputFcn`

• `OutputSelection`

• `Vectorization`

`matlab.ode.options.CVODESNonstiff Properties`

`matlab.ode.options.CVODESStiff Properties`

`matlab.ode.options.IDAS Properties`

• `InitialStep`

• `MaxStep`

• `MaxOrder`

Absolute error tolerance, specified as a positive scalar or vector. This tolerance is a threshold below which the value of the solution becomes unimportant. If the solution `|y|` is less than `AbsoluteTolerance`, then the solver does not need to obtain any correct digits in `|y|`. For this reason, keep in mind the scale of the solution components when setting the value of `AbsoluteTolerance`.

If `AbsoluteTolerance` is a vector, then it must be the same length as the number of solution components. If `AbsoluteTolerance` is a scalar, then the value applies to all solution components.

At each step, the ODE solver estimates the local error `e` in the `i`th component of the solution. To be successful, the step must have an acceptable error, as determined by both the relative and absolute error tolerances:

`|e(i)| <= max(RelativeTolerance*abs(y(i)),AbsoluteTolerance(i))`

Data Types: `single` | `double`

Relative error tolerance, specified as a positive scalar. This tolerance measures the error relative to the magnitude of each solution component. The relative error tolerance controls the number of correct digits in all solution components, except those smaller than `AbsoluteTolerance`.

At each step, the ODE solver estimates the local error `e` in the `i`th component of the solution. To be successful, the step must have an acceptable error, as determined by both the relative and absolute error tolerances:

`|e(i)| <= max(RelativeTolerance*abs(y(i)),AbsoluteTolerance(i))`

Data Types: `single` | `double`

## Object Functions

 `solve` Solve ODE over interval or at specified points `solutionFcn` Construct function that interpolates ODE solution `selectSolver` Change ODE solver

## Examples

collapse all

Create an empty `ode` object, and then specify values for the `ODEFcn` and `InitialValue` properties.

```F = ode; F.ODEFcn = @(t,y) 2*t; F.InitialValue = 0;```

Check which solver is selected for the problem, and then solve the equation over the time range `[0 10]`.

`F.SelectedSolver`
```ans = SolverID enumeration ode45 ```
`sol = solve(F,0,10)`
```sol = ODEResults with properties: Time: [0 0.2500 0.5000 0.7500 1 1.2500 1.5000 1.7500 2 2.2500 2.5000 2.7500 3 3.2500 3.5000 3.7500 4 4.2500 4.5000 4.7500 5 5.2500 5.5000 5.7500 6 6.2500 6.5000 6.7500 7 7.2500 7.5000 7.7500 8 8.2500 8.5000 8.7500 9 9.2500 9.5000 9.7500 10] Solution: [0 0.0625 0.2500 0.5625 1.0000 1.5625 2.2500 3.0625 4 5.0625 6.2500 7.5625 9 10.5625 12.2500 14.0625 16 18.0625 20.2500 22.5625 25 27.5625 30.2500 33.0625 36 39.0625 42.2500 45.5625 49 52.5625 56.2500 60.0625 64 68.0625 ... ] (1x41 double) ```

Plot the results.

`plot(sol.Time,sol.Solution,"-o")`

The Van der Pol oscillator equation is a second-order differential equation. The equation includes a parameter $\mu$, and the equation becomes stiff when the value of $\mu$ is large.

`$\frac{{d}^{2}x}{d{t}^{2}}-\mu \left(1-{x}^{2}\right)\frac{dx}{dt}+x=0$`

Using the substitutions ${\mathit{y}}_{1}=\mathit{x}$ and ${\mathit{y}}_{2}=\frac{\mathrm{dx}}{\mathrm{dt}}$ produces a system of two first-order equations.

`$\begin{array}{l}\frac{d{y}_{1}}{dt}={y}_{2}\\ \frac{d{y}_{2}}{dt}=\mu \left(1-{y}_{1}^{2}\right){y}_{2}-{y}_{1}\end{array}$`

The Jacobian matrix for these equations is the matrix of partial derivatives of each equation with respect to both ${\mathit{y}}_{1}$ and ${\mathit{y}}_{2}$.

`$J=\left[\begin{array}{cc}\frac{\partial {f}_{1}}{\partial {y}_{1}}& \frac{\partial {f}_{1}}{\partial {y}_{2}}\\ \frac{\partial {f}_{2}}{\partial {y}_{1}}& \frac{\partial {f}_{2}}{\partial {y}_{2}}\end{array}\right]=\left[\begin{array}{cc}0& 1\\ -2\mu {y}_{1}{y}_{2}-1& \mu \left(1-{y}_{1}^{2}\right)\end{array}\right]$`

Solve the Van der Pol oscillator using $\mu =1000$ and initial values of `[2; 0]` by creating an `ode` object to represent the problem.

• Store the value of $\mu$ in the `Parameters` property.

• Specify the initial values in the `InitialValue` property.

• Specify the system of equations in the `ODEFcn` property, specifying three input arguments so that the value for $\mu$ is passed to the function.

• Specify a function that calculates the Jacobian matrix in the `Jacobian` property, specifying three input arguments so that the value for $\mu$ is passed to the function.

```F = ode; F.Parameters = 1000; F.InitialValue = [2; 0]; F.ODEFcn = @(t,y,p) [y(2); p(1)*(1-y(1)^2)*y(2)-y(1)]; F.Jacobian = @(t,y,p) [0 1; -2*p(1)*y(1)*y(2)-1 p(1)*(1-y(1)^2)];```

Display the `ode` object. The `SelectedSolver` property shows that the `ode15s` solver was automatically chosen for this problem.

`F`
```F = ode with properties: Problem definition ODEFcn: @(t,y,p)[y(2);p(1)*(1-y(1)^2)*y(2)-y(1)] Parameters: 1000 InitialTime: 0 InitialValue: [2x1 double] Jacobian: [1x1 odeJacobian] EquationType: standard Solver properties AbsoluteTolerance: 1.0000e-06 RelativeTolerance: 1.0000e-03 Solver: auto SelectedSolver: ode15s Show all properties ```

Solve the system of equations over the time interval `[0 3000]` by using the `solve` method. Plot the first solution component.

```S = solve(F,0,3000); plot(S.Time,S.Solution(1,:),"-o")```

Create an `ode` object for the van der Pol equations with `mu = 1000` using the built-in function file `vdp1000.m`. Specify the initial values of `y'` and `y''` as `2` and `0`, respectively.

`F = ode(ODEFcn=@vdp1000,InitialValue=[2;0])`
```F = ode with properties: Problem definition ODEFcn: @vdp1000 InitialTime: 0 InitialValue: [2x1 double] EquationType: standard Solver properties AbsoluteTolerance: 1.0000e-06 RelativeTolerance: 1.0000e-03 Solver: auto SelectedSolver: ode45 Show all properties ```

The automatically selected solver is `ode45`. Solve the ODE system over the time interval `[0 3000]`. Measure the time the `ode45` solver takes by using `tic` and `toc`.

```tic sol45 = solve(F,0,3000); toc```
```Elapsed time is 11.110268 seconds. ```

Call the `selectSolver` function with the `DetectStiffness` name-value argument to choose a solver using a stiffness detection heuristic. Also specify the `IntervalLength` name-value argument as `3000`.

`F.Solver = selectSolver(F,DetectStiffness="on",IntervalLength=3000)`
```F = ode with properties: Problem definition ODEFcn: @vdp1000 InitialTime: 0 InitialValue: [2x1 double] Jacobian: [] EquationType: standard Solver properties AbsoluteTolerance: 1.0000e-06 RelativeTolerance: 1.0000e-03 Solver: ode15s Show all properties ```

The selected solver is now `ode15s`. Solving the ODE system over the time interval `[0 3000]` using the `ode15s` solver is faster and more efficient than using the `ode45` solver.

```tic sol15s = solve(F,0,3000); toc```
```Elapsed time is 0.243462 seconds. ```

Consider this system of first-order equations.

`$\begin{array}{l}{y}_{1}^{\prime }={y}_{1}{y}_{3}-{y}_{2}\\ {y}_{2}^{\prime }={y}_{1}-1\\ 0={y}_{1}+{y}_{2}+{y}_{3}\end{array}$`

The left side of the equations contain time derivatives, ${\mathit{y}}_{\mathit{n}}^{\prime }=\frac{{\mathrm{dy}}_{\mathit{n}}}{\mathrm{dt}}$. However, because the derivative for ${\mathit{y}}_{3}$ does not appear in the system, the equations define a system of differential algebraic equations. Rewriting the system in the form ${\mathrm{My}}^{\prime }=\mathit{f}\left(\mathit{t},\mathit{y}\right)$ shows a constant, singular mass matrix on the left side.

`$\left[\begin{array}{ccc}1& 0& 0\\ 0& 1& 0\\ 0& 0& 0\end{array}\right]\underset{}{\overset{˙}{\left[\begin{array}{c}{y}_{1}^{\prime }\\ {y}_{2}^{\prime }\\ {y}_{3}^{\prime }\end{array}\right]}}=\left[\begin{array}{c}{y}_{1}{y}_{3}-{y}_{2}\\ {y}_{1}-1\\ {y}_{1}+{y}_{2}+{y}_{3}\end{array}\right]$`

Solve the system of equations using the initial values `[1 1 -2]` by creating an `ode` object to represent the problem.

• Specify the initial values in the `InitialValue` property.

• Specify the system of equations as an anonymous function in the `ODEFcn` property.

• Use an `odeMassMatrix` object to specify the constant, singular mass matrix in the `MassMatrix` property.

```F = ode; F.InitialValue = [1 1 -2]; F.ODEFcn = @(t,y) [y(1)*y(3)-y(2); y(1)-1; y(1)+y(2)+y(3)]; F.MassMatrix = odeMassMatrix(MassMatrix=[1 0 0; 0 1 0; 0 0 0],Singular="yes");```

Display the ode object. The `SelectedSolver` property shows that the `ode15s` solver was automatically chosen for this problem.

`F`
```F = ode with properties: Problem definition ODEFcn: @(t,y)[y(1)*y(3)-y(2);y(1)-1;y(1)+y(2)+y(3)] InitialTime: 0 InitialValue: [1 1 -2] Jacobian: [] MassMatrix: [1x1 odeMassMatrix] EquationType: standard Solver properties AbsoluteTolerance: 1.0000e-06 RelativeTolerance: 1.0000e-03 Solver: auto SelectedSolver: ode15s Show all properties ```

Solve the system of equations over the time interval `[0 10]` by using the `solve` method. Plot all three solution components.

```S = solve(F,0,10); plot(S.Time,S.Solution,"-o") legend("y_1","y_2","y_3",Location="southeast")```

Solve an ODE system with two equations and two parameters, and perform sensitivity analysis on the parameters.

Create an `ode` object to represent this system of equations.

`$\begin{array}{l}\frac{{\mathrm{dy}}_{1}}{\mathrm{dt}}={\mathit{p}}_{1}{\mathit{y}}_{1}-{\mathit{y}}_{2}\\ \frac{{\mathrm{dy}}_{2}}{\mathrm{dt}}={-\mathit{p}}_{2}{\mathit{y}}_{2}\end{array}$`

Specify the initial conditions as ${\mathit{y}}_{1}\left(0\right)=2$ and ${\mathit{y}}_{2}\left(0\right)=3$, and parameter values of ${\mathit{p}}_{1}=0.05$ and ${\mathit{p}}_{2}=1.5$. To enable sensitivity analysis of the parameters, set the `Sensitivity` property of the `ode` object to an `odeSensitivity` object.

```p = [0.05 1.5]; F = ode(ODEFcn=@(t,y,p) [p(1)*y(1)-y(2); -p(2)*y(2)], ... InitialValue=[2 3], ... Parameters=p, ... Sensitivity=odeSensitivity)```
```F = ode with properties: Problem definition ODEFcn: @(t,y,p)[p(1)*y(1)-y(2);-p(2)*y(2)] Parameters: [0.0500 1.5000] InitialTime: 0 InitialValue: [2 3] Sensitivity: [1x1 odeSensitivity] EquationType: standard Solver properties AbsoluteTolerance: 1.0000e-06 RelativeTolerance: 1.0000e-03 Solver: auto SelectedSolver: cvodesnonstiff Show all properties ```

Because the equations are nonstiff and sensitivity analysis is enabled, the `ode` object automatically chooses the `cvodesnonstiff` solver for this problem.

Solve the ODE over the time interval `[0 5]`, and plot the solution for each component.

`S = solve(F,0,5)`
```S = ODEResults with properties: Time: [0 2.9540e-09 2.9543e-05 2.2465e-04 4.1976e-04 0.0024 0.0080 0.0137 0.0245 0.0353 0.0611 0.0869 0.1499 0.2129 0.3169 0.4208 0.5248 0.7204 0.9161 1.1118 1.3075 1.5031 1.6988 1.8945 2.0901 2.2858 2.4815 2.6772 2.8728 ... ] (1x40 double) Solution: [2x40 double] Sensitivity: [2x2x40 double] ```
```plot(S.Time,S.Solution(1,:),"-o",S.Time,S.Solution(2,:),"-o") legend("y1","y2")```

The values in `S.Sensitivity` are partial derivatives of the equations with respect to the parameters. To examine the effects of the parameter values during the integration, plot the sensitivity values.

```figure hold on plot(S.Time,squeeze(S.Sensitivity(1,1,:)),"-o") plot(S.Time,squeeze(S.Sensitivity(1,2,:)),"-o") plot(S.Time,squeeze(S.Sensitivity(2,1,:)),"-o") plot(S.Time,squeeze(S.Sensitivity(2,2,:)),"-o") legend("p1,eq1","p2,eq1","p1,eq2","p2,eq2") hold off```

Consider a ball thrown upward with an initial velocity $\mathrm{dy}/\mathrm{dt}$. The ball is subject to acceleration due to gravity aimed downward, so its acceleration is

`$\frac{{d}^{2}y}{d{t}^{2}}=-g.$`

Rewriting the equation as a first-order system of equations with the substitutions ${\mathit{y}}_{1}=\mathit{y}$ and ${\mathit{y}}_{2}=\frac{\mathrm{dy}}{\mathrm{dt}}$ yields

`$\begin{array}{l}{y}_{1}^{\prime }={y}_{2}\\ {y}_{2}^{\prime }=-g.\end{array}$`

Solve the equations for the position ${\mathit{y}}_{1}$ and velocity ${\mathit{y}}_{2}$ of the ball over time.

Define Equations and Initial Conditions

Create a function handle for the first-order system of equations that accepts two inputs for `(t,y)`. Use the value $\mathit{g}=9.8$ for the acceleration due to gravity.

`dydt = @(t,y) [y(2); -9.8];`

Next, create a vector with the initial conditions. The ball starts at position ${\mathit{y}}_{1}=3$ at $\mathit{t}=0$ as it is thrown upward with initial velocity ${\mathit{y}}_{2}=20$.

`y0 = [3 20];`

Model Ball Bounces as Events

The ball initially travels upward until the force due to gravity causes it to change direction and head back down to the ground. If you solve the equations without more consideration, then the ball falls back downward forever without striking the ground. Instead, you can use an event function to detect when the position of the ball goes to zero where the ground is located. Because the solution component ${\mathit{y}}_{1}=\mathit{y}$ is the position of the ball, the event function tracks the value of ${\mathit{y}}_{1}$ so that an event triggers whenever ${\mathit{y}}_{1}=0$.

Create a function handle for the event function that accepts two inputs for `(t,y)`.

`bounceEvent = @(t,y) y(1);`

When the ball strikes the ground, its direction changes again as it heads back upwards with a new (smaller) initial velocity. To model this situation, use a callback function along with the event function. When an event triggers, the ODE solver invokes the callback function. The callback function resets the position and initial velocity of the ball so that the integration can restart with the correct initial conditions. When an event occurs, the callback function sets the position ${\mathit{y}}_{1}=0$ and attenuates the velocity by a factor of $0.9$ while reversing its direction back upward. Define a callback function that performs these actions.

```function [stop,y] = bounceResponse(t,y) stop = false; y(1) = 0; y(2) = -0.9*y(2); end ```

(The callback function is included as a local function at the end of the example.)

Create an `odeEvent` object to represent the bouncing ball events. Specify `Direction="descending"` so that only events where the position is decreasing are detected. Also, specify `Response="callback"` so that the solver invokes the callback function when an event occurs.

```E = odeEvent(EventFcn=bounceEvent, ... Direction="descending", ... Response="callback", ... CallbackFcn=@bounceResponse)```
```E = odeEvent with properties: EventFcn: @(t,y)y(1) Direction: descending Response: callback CallbackFcn: @bounceResponse ```

Solve Equations

Create an `ode` object for the problem, specifying the equations `dydt`, initial conditions `y0`, and events `E` as property values.

`F = ode(ODEFcn=dydt,InitialValue=y0,EventDefinition=E);`

Integrate the equations over the time interval `[0 30]` by using the `solve` method. Specify `Refine=8` to generate 8 points per step. The resulting object has properties for the time and solution, and because events are being tracked, the object also displays properties related to the events that triggered during the integration.

`S = solve(F,0,30,Refine=8)`
```S = ODEResults with properties: Time: [0 0.0038 0.0075 0.0113 0.0151 0.0188 0.0226 0.0264 0.0301 0.0490 0.0678 0.0867 0.1055 0.1243 0.1432 0.1620 0.1809 0.2751 0.3692 0.4634 0.5576 0.6518 0.7460 0.8402 0.9344 1.3094 1.6844 2.0594 2.4344 2.8094 3.1844 ... ] (1x468 double) Solution: [2x468 double] EventTime: [4.2265 8.1607 11.7015 14.8882 17.7563 20.3375 22.6606 24.7514 26.6331 28.3267 29.8509] EventSolution: [2x11 double] EventIndex: [1 1 1 1 1 1 1 1 1 1 1] ```

Plot Results

Plot the position ${\mathit{y}}_{1}$ of the ball over time, marking the initial position with a green circle and events with red circles.

```plot(S.Time,S.Solution(1,:),"--") hold on plot(S.EventTime,S.EventSolution(1,:),"ro") plot(0,y0(1),"go") hold off ylim([0 25]) xlabel("Time") ylabel("Position y_1")```

Local Functions

```function [stop,y] = bounceResponse(t,y) stop = false; y(1) = 0; y(2) = -0.9*y(2); end```

## Version History

Introduced in R2023b

expand all