Documentation

# ode23t

Solve moderately stiff ODEs and DAEs — trapezoidal rule

## Syntax

``````[t,y] = ode23t(odefun,tspan,y0)``````
``````[t,y] = ode23t(odefun,tspan,y0,options)``````
``````[t,y,te,ye,ie] = ode23t(odefun,tspan,y0,options)``````
``sol = ode23t(___)``

## Description

example

``````[t,y] = ode23t(odefun,tspan,y0)```, where `tspan = [t0 tf]`, integrates the system of differential equations $y\text{'}=f\left(t,y\right)$ from `t0` to `tf` with initial conditions `y0`. Each row in the solution array `y` corresponds to a value returned in column vector `t`.All MATLAB® ODE solvers can solve systems of equations of the form $y\text{'}=f\left(t,y\right)$, or problems that involve a mass matrix, $M\left(t,y\right)y\text{'}=f\left(t,y\right)$. The solvers all use similar syntaxes. The `ode23s` solver only can solve problems with a mass matrix if the mass matrix is constant. `ode15s` and `ode23t` can solve problems with a mass matrix that is singular, known as differential-algebraic equations (DAEs). Specify the mass matrix using the `Mass` option of `odeset`.```

example

``````[t,y] = ode23t(odefun,tspan,y0,options)``` also uses the integration settings defined by `options`, which is an argument created using the `odeset` function. For example, use the `AbsTol` and `RelTol` options to specify absolute and relative error tolerances, or the `Mass` option to provide a mass matrix.```
``````[t,y,te,ye,ie] = ode23t(odefun,tspan,y0,options)``` additionally finds where functions of (t,y), called event functions, are zero. In the output, `te` is the time of the event, `ye` is the solution at the time of the event, and `ie` is the index of the triggered event.For each event function, specify whether the integration is to terminate at a zero and whether the direction of the zero crossing matters. Do this by setting the `'Events'` property to a function, such as `myEventFcn` or `@myEventFcn`, and creating a corresponding function: [`value`,`isterminal`,`direction`] = `myEventFcn`(`t`,`y`). For more information, see ODE Event Location.```
````sol = ode23t(___)` returns a structure that you can use with `deval` to evaluate the solution at any point on the interval `[t0 tf]`. You can use any of the input argument combinations in previous syntaxes.```

## Examples

collapse all

Simple ODEs that have a single solution component can be specified as an anonymous function in the call to the solver. The anonymous function must accept two inputs `(t,y)` even if one of the inputs is not used.

Solve the ODE

`${y}^{\prime }=-10t.$`

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

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

Plot the solution.

`plot(t,y,'-o')`

An example of a stiff system of equations is the van der Pol equations in relaxation oscillation. The limit cycle has regions where the solution components change slowly and the problem is quite stiff, alternating with regions of very sharp change where it is not stiff.

The system of equations is:

The initial conditions are and . The function `vdp1000` ships with MATLAB® and encodes the equations.

```function dydt = vdp1000(t,y) %VDP1000 Evaluate the van der Pol ODEs for mu = 1000. % % See also ODE15S, ODE23S, ODE23T, ODE23TB. % Jacek Kierzenka and Lawrence F. Shampine % Copyright 1984-2014 The MathWorks, Inc. dydt = [y(2); 1000*(1-y(1)^2)*y(2)-y(1)]; ```

Solving this system using `ode45` with the default relative and absolute error tolerances (`1e-3` and `1e-6`, respectively) is extremely slow, requiring several minutes to solve and plot the solution. `ode45` requires millions of time steps to complete the integration, due to the areas of stiffness where it struggles to meet the tolerances.

This is a plot of the solution obtained by `ode45`, which takes a long time to compute. Notice the enormous number of time steps required to pass through areas of stiffness.

Solve the stiff system using the `ode23t` solver, and then plot the first column of the solution `y` against the time points `t`. The `ode23t` solver passes through stiff areas with far fewer steps than `ode45`.

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

`ode23t` only works with functions that use two input arguments, `t` and `y`. However, you can pass in extra parameters by defining them outside the function and passing them in when you specify the function handle.

Solve the ODE

Rewriting the equation as a first-order system yields

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

```function dydt = odefcn(t,y,A,B) dydt = zeros(2,1); dydt(1) = y(2); dydt(2) = (A/B)*t.*y(1); ```

Solve the ODE using `ode23t`. Specify the function handle such that it passes in the predefined values for `A` and `B` to `odefcn`.

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

Plot the results.

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

The `ode15s` solver is a good first choice for most stiff problems. However, the other stiff solvers might be more efficient for certain types of problems. This example solves a stiff test equation using all four stiff ODE solvers.

Consider the test equation

`${y}^{\prime }=-\lambda y.$`

The equation becomes increasingly stiff as the magnitude of $\lambda$ increases. Use $\lambda =1×1{0}^{9}$ and the initial condition $y\left(0\right)=1$ over the time interval `[0 0.5]`. With these values, the problem is stiff enough that `ode45` and `ode23` struggle to integrate the equation. Also, use `odeset` to pass in the constant Jacobian $J=\frac{\partial f}{\partial y}=-\lambda$ and turn on the display of solver statistics.

```lambda = 1e9; y0 = 1; tspan = [0 0.5]; opts = odeset('Jacobian',-lambda,'Stats','on');```

Solve the equation with `ode15s`, `ode23s`, `ode23t`, and `ode23tb`. Make subplots for comparison.

```subplot(2,2,1) tic, ode15s(@(t,y) -lambda*y, tspan, y0, opts), toc```
```104 successful steps 1 failed attempts 212 function evaluations 0 partial derivatives 21 LU decompositions 210 solutions of linear systems Elapsed time is 1.951681 seconds. ```
```title('ode15s') subplot(2,2,2) tic, ode23s(@(t,y) -lambda*y, tspan, y0, opts), toc```
```63 successful steps 0 failed attempts 191 function evaluations 0 partial derivatives 63 LU decompositions 189 solutions of linear systems Elapsed time is 0.484418 seconds. ```
```title('ode23s') subplot(2,2,3) tic, ode23t(@(t,y) -lambda*y, tspan, y0, opts), toc```
```95 successful steps 0 failed attempts 125 function evaluations 0 partial derivatives 28 LU decompositions 123 solutions of linear systems Elapsed time is 0.961416 seconds. ```
```title('ode23t') subplot(2,2,4) tic, ode23tb(@(t,y) -lambda*y, tspan, y0, opts), toc```
```71 successful steps 0 failed attempts 167 function evaluations 0 partial derivatives 23 LU decompositions 236 solutions of linear systems Elapsed time is 0.578629 seconds. ```
`title('ode23tb')`

The stiff solvers all perform well, but `ode23s` completes the integration with the fewest steps and runs the fastest for this particular problem. Since the constant Jacobian is specified, none of the solvers need to calculate partial derivatives to compute the solution. Specifying the Jacobian benefits `ode23s` the most since it normally evaluates the Jacobian in every step.

For general stiff problems, the performance of the stiff solvers varies depending on the format of the problem and specified options. Providing the Jacobian matrix or sparsity pattern always improves solver efficiency for stiff problems. But since the stiff solvers use the Jacobian differently, the improvement can vary significantly. Practically speaking, if a system of equations is very large or needs to be solved many times, then it is worthwhile to investigate the performance of the different solvers to minimize execution time.

This example shows how to use `ode23t` to solve a stiff differential algebraic equation (DAE) problem that arises from an electrical circuit. This problem is originally from p.377 of E. Hairer and G. Wanner, Solving Ordinary Differential Equations II: Stiff and Differential-Algebraic Problems, 2nd ed. (Berlin: Springer, 1996).

The one-transistor amplifier problem coded by amp1dae.m can be rewritten in semi-explicit form, but this example solves it in its original form . The problem includes a constant, singular mass matrix that is not diagonal. The transistor amplifier circuit contains six resistors, three capacitors, and a transistor.

The initial voltage signal is , and the other parameters are constant. The goal is to solve for the output voltage through node 5.

Set the values of the initial and operating voltages.

```Ue = @(t) 0.4*sin(200*pi*t); Ub = 6; ```

Using Kirchoff's law to equalize the current through each node (1 through 5), you obtain a system of five equations. The mass matrix of this system has the form

where for .

Use the `odeset` function to pass in the mass matrix. Leave the `'MassSingular'` option at its default value of `'maybe'` to test the automatic detection of a DAE problem by the solver.

```c = 1e-6 * (1:3); M = zeros(5,5); M(1,1) = -c(1); M(1,2) = c(1); M(2,1) = c(1); M(2,2) = -c(1); M(3,3) = -c(2); M(4,4) = -c(3); M(4,5) = c(3); M(5,4) = c(3); M(5,5) = -c(3); opts = odeset('Mass',M); ```

The function `transampdae.m` contains the system of equations for this example. Save this function in your current folder to run the example.

```function dudt = transampdae(t,u) % Define voltages and constant parameters Ue = @(t) 0.4*sin(200*pi*t); Ub = 6; R0 = 1000; R15 = 9000; alpha = 0.99; beta = 1e-6; Uf = 0.026; f23 = beta*(exp((u(2) - u(3))/Uf) - 1); dudt = [ -(Ue(t) - u(1))/R0 -(Ub/R15 - u(2)*2/R15 - (1-alpha)*f23) -(f23 - u(3)/R15) -((Ub - u(4))/R15 - alpha*f23) (u(5)/R15) ]; end ```

Set the initial conditions. This example uses the consistent initial conditions computed by Hairer and Wanner.

```u0(1) = 0; u0(2) = Ub/2; u0(3) = Ub/2; u0(4) = Ub; u0(5) = 0; ```

Solve the DAE system over the time interval `[0 0.05]` using `ode23t`.

```tspan = [0 0.05]; [t,u] = ode23t(@transampdae,tspan,u0,opts); ```

Plot the initial voltage and output voltage .

```plot(t,Ue(t),'o',t,u(:,5),'.') axis([0 0.05 -3 2]); legend('Input Voltage U_e(t)','Output Voltage U_5(t)', 'Location', 'NorthWest'); title('One transistor amplifier DAE problem solved by ODE23T'); xlabel('t'); ```

## Input Arguments

collapse all

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

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

For example, to solve $y\text{'}=5y-3$, use the function:

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

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

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

use the function:

```function dydt = odefun(t,y) dydt = zeros(2,1); dydt(1) = y(1)+2*y(2); dydt(2) = 3*y(1)+2*y(2);```

For information on how to provide additional parameters to the function `odefun`, see Parameterizing Functions.

Example: `@myFcn`

Data Types: `function_handle`

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

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

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

• If `tspan` has more than two elements `[t0,t1,t2,...,tf]`, then the solver returns the solution evaluated at the given points. However, the solver does not step precisely to each point specified in `tspan`. Instead, the solver uses its own internal steps to compute the solution, then evaluates the solution at the requested points in `tspan`. The solutions produced at the specified points are of the same order of accuracy as the solutions computed at each internal step.

Specifying several intermediate points has little effect on the efficiency of computation, but for large systems it can affect memory management.

The values of `tspan` are used by the solver to calculate suitable values for `InitialStep` and `MaxStep`:

• If `tspan` contains several intermediate points `[t0,t1,t2,...,tf]`, then the specified points give an indication of the scale for the problem, which can affect the value of `InitialStep` used by the solver. Therefore, the solution obtained by the solver might be different depending on whether you specify `tspan` as a two-element vector or as a vector with intermediate points.

• The initial and final values in `tspan` are used to calculate the maximum step size `MaxStep`. Therefore, changing the initial or final values in `tspan` could lead to the solver using a different step sequence, which might change the solution.

Example: `[1 10]`

Example: ```[1 3 5 7 9 10]```

Data Types: `single` | `double`

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

Data Types: `single` | `double`

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

Example: `options = odeset('RelTol',1e-5,'Stats','on','OutputFcn',@odeplot)` specifies a relative error tolerance of `1e-5`, turns on the display of solver statistics, and specifies the output function `@odeplot` to plot the solution as it is computed.

Data Types: `struct`

## Output Arguments

collapse all

Evaluation points, returned as a column vector.

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

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

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

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

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

Index of triggered event function, returned as a column vector. The event times in `te` correspond to the solutions returned in `ye`, and `ie` specifies which event occurred.

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

Structure FieldDescription

`sol.x`

Row vector of the steps chosen by the solver.

`sol.y`

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

`sol.solver`

Solver name.

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

Structure FieldDescription

`sol.xe`

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

`sol.ye`

Solutions that correspond to events in `sol.xe`.

`sol.ie`

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

## Algorithms

`ode23t` is an implementation of the trapezoidal rule using a “free” interpolant. This solver is preferred over `ode15s` if the problem is only moderately stiff and you need a solution without numerical damping. `ode23t` also can solve differential algebraic equations (DAEs) [1], [2].

## References

[1] Shampine, L. F., M. W. Reichelt, and J.A. Kierzenka, “Solving Index-1 DAEs in MATLAB and Simulink”, SIAM Review, Vol. 41, 1999, pp. 538–552.

[2] Shampine, L. F. and M. E. Hosea, “Analysis and Implementation of TR-BDF2,” Applied Numerical Mathematics 20, 1996.