# numeric::odesolve

Numerical solution of an ordinary differential equation

### Use only in the MuPAD Notebook Interface.

This functionality does not run in MATLAB.

## Syntax

```numeric::odesolve(`f`, `t0 .. t`, `Y0`, <`method`>, <`RelativeError = rtol`>, <`AbsoluteError = atol`>, <`Stepsize = h`>, <`MaxStepsize = hmax`>, <`Alldata = n`>, <Symbolic>)
numeric::odesolve(`t0 .. t`, `f`, `Y0`, <`method`>, <`RelativeError = rtol`>, <`AbsoluteError = atol`>, <`Stepsize = h`>, <`MaxStepsize = hmax`>, <`Alldata = n`>, <Symbolic>)
```

## Description

```numeric::odesolve(f, t0..t, Y0)``` returns a numerical approximation of the solution Y(t) of the first order differential equation (dynamical system) , Y(t0) = Y0 with and .

`numeric::odesolve` is a general purpose solver able to deal with initial value problems of various kinds of ordinary differential equations. Equations of order p can be solved by `numeric::odesolve` after reformulation to dynamical system form. This can always be achieved by writing the equation as a first order system

for the vector . See Example 4.

The following single-step Runge-Kutta-type methods are implemented:

• `EULER1` (order 1)

• `RKF43` (order 3)

• `xRKF43` (order 3)

• `RKF34` (order 4)

• `xRKF34` (order 4)

• `RK4` (order 4)

• `RKF54a` (order 4)

• `RKF54b` (order 4)

• `DOPRI54` (order 4)

• `xDOPRI54` (order 4)

• `CK54` (order 4)

• `xRKF54a` (order 4)

• `xRKF54b` (order 4)

• `xCK54` (order 4)

• `RKF45a` (order 5)

• `RKF45b` (order 5)

• `DOPRI45` (order 5)

• `CK45` (order 5)

• `xRKF45a` (order 5)

• `xRKF45b` (order 5)

• `xDOPRI45` (order 5)

• `xCK45` (order 5)

• `DOPRI65` (order 5)

• `xDOPRI65` (order 5)

• `DOPRI56` (order 6)

• `xDOPRI56` (order 6)

• `BUTCHER6` (order 6)

• `RKF87` (order 7)

• `xRKF87` (order 7)

• `DOPRI87` (order 7)

• `xDOPRI87` (order 7)

• `RKF78` (order 8)

• `xRKF78` (order 8)

• `DOPRI78` (order 8)

• `xDOPRI78` (order 8)

• `GAUSS(s)` (order 2 s)

• `GAUSS` = `s`

For the Gauss methods, `GAUSS(s)` is equivalent to `GAUSS = s`. The positive integer `s` indicates the number of stages. The order of the `s` stage Gauss method is 2 s.

The utility function `numeric::ode2vectorfield` may be used to produce the input parameters f, t0, Y0 from a set of differential expressions representing the ODE. See Example 1.

The input data t0, t and Y0 must not contain symbolic objects which cannot be converted to floating point values via `float`. Numerical expressions such as , etc. are accepted.

The vector field f defining the dynamical system must be represented by a procedure with two input parameters: the scalar time t and the vector Y. `numeric::odesolve` internally calls this function with real floating-point values t and a list Y of floating-point values. It has to return the vector f(t, Y) either as a list or as a 1-dimensional array. The output of f may contain numerical expressions such as π, etc. However, all values must be convertible to real or complex floating point numbers by `float`.

Autonomous systems, where f(t, Y) does not depend on t, must also be represented by a procedure with 2 arguments `t` and `Y`.

Scalar functions `Y` also must be represented by a list or an array with one element. For instance, the input data for the scalar initial value problem may be of the form

 `f := proc(t,Y)` /* Y is a 1-dimensional vector */ `local y;` /* represented by a list with */ `begin` /* one element: Y = [y]. */ `   y := Y[1];` `   [t*sin(y)]` /* the output is a list with 1 element */ `end_proc:` `Y0 := [1]:` /* the initial value */

The numerical precision is controlled by the global variable `DIGITS`: an adaptive control of the step size keeps local relative discretization errors below `rtol=10^-DIGITS`, unless a different tolerance is specified via the option `RelativeError = rtol`. For small values of the solution vector Y, the absolute discretization error can be bounded by the threshold `atol` specified via the option `AbsoluteError = atol`.

If `AbsoluteError` is not specified, only relative discretization errors are controlled and kept below `rtol`.

The error control may be switched off by specifying a fixed ```Stepsize = h```.

 Note:   Only local errors are controlled by the adaptive mechanism. No control of the global error is provided!

With `Y := t -> numeric::odesolve(f, t_0..t, Y_0)`, the numerical solution can be repesented by a MuPAD® function: the call `Y(t)` will start the numerical integration from `t0` to `t`. A more sophisticated form of this function may be generated via ```Y := numeric::odesolve2(f, t0, Y0)```.

This equips `Y` with a remember mechanism that uses previously computed values to speed up the computation. See Example 2.

For systems of the special form with a matrix valued function f(t, Y), there is a special solver `numeric::odesolveGeometric` which preserves geometric features of the system more faithfully than `numeric::odesolve`.

## Environment Interactions

The function is sensitive to the environment variable `DIGITS`, which determines the numerical working precision.

## Examples

### Example 1

We compute the numerical solution y(10) of the initial value problem , y(0) = 2:

```f := proc(t, Y) begin [t*sin(Y[1])] end_proc: numeric::odesolve(f, 0..10, [2])```

Alternatively, the utility function `numeric::ode2vectorfield` can be used to generate the input parameters in a more intuitive way:

```[f, t0, Y0] := [numeric::ode2vectorfield({y'(t) = t*sin(y(t)), y(0) = 2}, [y(t)])]```

`numeric::odesolve(f, t0..10, Y0)`

`delete f, t0, Y0:`

### Example 2

We consider , y(0) = 1. The numerical solution may be represented by the function

`Y := t -> numeric::odesolve((t,Y) -> Y, 0..t, [1]):`

Calling `Y(t)` starts the numerical integration:

`Y(-5), Y(0), Y(1), Y(PI)`

`delete Y:`

### Example 3

We compute the numerical solution Y(π) = [x(π), y(π)] of the system

.

```f := (t, Y) -> [Y[1] + Y[2], Y[1] - Y[2]]: Y0 := [1, I]: numeric::odesolve(f, 0..PI, Y0)```

The solution of a linear dynamical system with a constant matrix A is . The solution of the system above can also be computed by:

```t := PI: tA := array(1..2, 1..2, [[t, t], [t, -t]]): numeric::expMatrix(tA, Y0)```

`delete f, Y0, t, tA:`

### Example 4

We compute the numerical solution y(1) of the differential equation with initial conditions y(0) = 0, . The second order equation is converted to a first order system for the vector :

.

```f := proc(t, Y) begin [Y[2], Y[1]^2] end_proc: Y0 := [0, 1]: numeric::odesolve(f, 0..1, Y0)```

`delete f, Y0:`

### Example 5

We demonstrate how numerical data can be obtained on a user defined time mesh `t[i]`. The initial value problem is , y(0) = 1, the sample points are t0 = 0.0, t1 = 0.1, …, t100 = 10.0. First, we define the differential equation and the initial condition:

```f := (t, Y) -> [sin(t) - Y[1]]: Y[0] := [1]:```

We define the time mesh:

`for i from 0 to 100 do t[i] := i/10 end_for:`

The equation is integrated iteratively from `t[i-1]` to `t[i]` with a working precision of 4 significant decimal places:

```DIGITS := 4: for i from 1 to 100 do Y[i] := numeric::odesolve(f, t[i-1]..t[i], Y[i-1]) end_for:```

The following mesh data are produced:

`[t[i], Y[i]] \$ i = 0..100`

These data can be displayed by a list plot:

```plotpoints := [[t[i], op(Y[i])] \$ i = 0..100]: plot(plot::PointList2d(plotpoints, PointColor = RGB::Black)):```

The same plot can be obtained directly via `plot::Ode2d`:

```plot(plot::Ode2d( [t[i] \$ i = 0..100], f, Y[0], [(t, Y) -> [t, Y[1]], Style = Points, Color = RGB::Black]))```

`delete f, t, DIGITS, Y, plotpoints:`

### Example 6

We compute the numerical solution y(1) of by the classical 4-th order Runge-Kutta method `RK4`. By internal local extrapolation, its effective order is 5:

```f := (t, Y) -> Y: DIGITS := 13: numeric::odesolve(f, 0..1, [1], RK4)```

Next, we use local extrapolation `xRKF78` of the 8-th order submethod of the Runge-Kutta-Fehlberg pair `RKF78`. This scheme has effective order 9:

`numeric::odesolve(f, 0..1, [1], xRKF78)`

Both methods yield the same result because of the internal adaptive error control. However, due to its higher order, the method `xRKF78` is faster.

`delete f, DIGITS:`

### Example 7

We consider the stiff ODE . The default method `DOPRI78` is explicit and not very efficient in solving very stiff problems:

```f := (t, Y) -> [10^4*(cos(t) - Y[1])]: t0 := time(): numeric::odesolve(f, 0..1, [1]), (time() - t0)*msec```

We use the implicit A-stable method `GAUSS(6)`. For this stiff problem, it is more efficient than the default method `DOPRI78`:

```t0 := time(): numeric::odesolve(f, 0..1, [1], GAUSS(6)), (time() - t0)*msec```

`delete t0:`

### Example 8

We consider the initial value problem , y(0) = 1. We note that the numerical evaluation of the right hand side of the equation suffers from cancellation effects, when |y| is small.

```f := (t, Y) -> [-10^20*Y[1]*(1 - cos(Y[1]))]: Y0 := [1]:```

We first attempt to compute y(1) with a working precision of 6 digits using the default setting ```RelativeError = 10^-DIGITS=10^(-6)```:

`DIGITS := 6: numeric::odesolve(f, 0..1, Y0)`

Due to numerical round-off in the internal steps, no digit of this result is correct. Next, we use a working precision of 20 significant decimal places to eliminate roundoff effects:

```DIGITS := 20: numeric::odesolve(f, 0..1, Y0, RelativeError = 10^(-6))```

Since relative local discretization errors are of the magnitude 10- 6, not all displayed digits are trustworthy. We finally use a working precision of 20 digits and constrain the local relative discretization errors by the tolerance 10- 10:

`numeric::odesolve(f, 0..1, Y0, RelativeError = 10^(-10))`

`delete f, Y0, DIGITS:`

### Example 9

We compute the numerical solution y(1) of , y(0) = 1 with various methods and various constant step sizes. We compare the result with the exact solution .

```f := (t, Y) -> Y: Y0 := [1]:```

We first use the Euler method of order 1 with two different step sizes:

```Y := numeric::odesolve(f, 0..1, Y0, EULER1, Stepsize = 0.1): Y, globalerror = float(exp(1)) - Y[1]```

Decreasing the step size by a factor of 10 should reduce the global error by a factor of 10. Indeed:

```Y := numeric::odesolve(f, 0..1, Y0, EULER1, Stepsize = 0.01): Y, globalerror = float(exp(1)) - Y[1]```

Next, we use the classical Runge-Kutta method of order 4 with two different step sizes:

```Y := numeric::odesolve(f, 0..1, Y0, RK4, Stepsize = 0.1): Y, globalerror = float(exp(1)) - Y[1]```

Decreasing the step size by a factor of 10 in a 4-th order scheme should reduce the global error by a factor of 104. Indeed:

```Y := numeric::odesolve(f, 0..1, Y0, RK4, Stepsize = 0.01): Y, globalerror = float(exp(1)) - Y[1]```

`delete f, Y0, Y:`

### Example 10

We integrate , y(0) = 1 over the interval t ∈ [0, 0.99] with a working precision of 4 digits. The exact solution is . Note the singularity at t = 1.

```DIGITS := 4: f := (t, Y) -> [Y[1]^2]: Y0 := [1]:```

The option `Alldata`, equivalent to ```Alldata = 1```, yields all mesh data generated during the internal adaptive process:

`numeric::odesolve(f, 0..0.99, Y0, Alldata)`

With `Alldata = 2`, only each second point is returned:

`numeric::odesolve(f, 0..0.99, Y0, Alldata = 2)`

One can control the time mesh using the option ```Stepsize = h```:

`numeric::odesolve(f, 0..0.99, Y0, Stepsize=0.1, Alldata = 1)`

However, with the option `Stepsize = h`, no automatic error control is provided by `numeric::odesolve`. Note the poor approximation y(t) = 94.3 for t = 0.99 (the exact value is y(0.99) = 100.0). The next computation with smaller step size yields better results:

`numeric::odesolve(f, 0..0.99, Y0, Stepsize = 0.01, Alldata = 10)`

Example 5 demonstrates how accurate numerical data on a user defined time mesh can be generated using the automatic error control of `numeric::odesolve`.

`delete DIGITS, f, Y0:`

### Example 11

The second order equation is written as the dynamical system , for the vector Y = [y, z]. A single symbolic step

of the Euler method is computed:

```f := proc(t, Y) begin [Y[2], -sin(Y[1])] end_proc: numeric::odesolve(f, t0..t0+h, [y0, z0], EULER1, Symbolic)```

`delete f:`

## Parameters

 `f` A procedure representing the vector field `t0` A numerical real value for the initial time `t` A numerical real value (the "time") `Y0` A list or 1-dimensional array of numerical values representing the initial condition `method` One of the Runge-Kutta schemes listed below.

## Options

`BUTCHER6`, `CK45`, `CK54`, `DOPRI45`, `DOPRI54`, `DOPRI56`, `DOPRI65`, `DOPRI78`, `DOPRI87`, `EULER1`, `RK4`, `RKF34`, `RKF43`, `RKF45a`, `RKF45b`, `RKF54a`, `RKF54b`, `RKF78`, `RKF87`, `xCK45`, `xCK54`, `xDOPRI45`, `xDOPRI54`, `xDOPRI56`, `xDOPRI65`, `xDOPRI78`, `xDOPRI87`, `xRKF34`, `xRKF43`, `xRKF45a`, `xRKF45b`, `xRKF54a`, `xRKF54b`, `xRKF78`, `xRKF87`

Name of the Runge-Kutta scheme. See Example 6. For details on these schemes, see the Algorithms section.

`GAUSS`

Name of the Runge-Kutta scheme specified as `GAUSS(s)` or ```GAUSS = s```.

The methods `GAUSS(s)` or, equivalently, ```GAUSS = s``` are the implicit Gauss methods with `s` stages of order 2 s.

These methods are implicit A-stable schemes. The time steps are rather costly to compute. The Gauss methods are useful for integrating stiff ODEs. For non-stiff ODEs, there is usually no need to change the default method `DOPRI78`. This method is an embedded Runge-Kutta pair of orders 7 and 8.

Further, the Gauss methods are symplectic methods. When used with constant step size (`Stepsize = h`), numerical integration of Hamiltonian systems benefits from this property.

See Example 7.

`RelativeError`, `AbsoluteError`

Option specified as `RelativeError = rtol` forces internal numerical Runge-Kutta steps to use step sizes with relative local discretization errors below `rtol`. This tolerance must be a positive numerical real value not smaller than . The default tolerance is `RelativeError = 10^(-DIGITS)`.

Option specified as `AbsoluteError = atol` forces internal numerical Runge-Kutta steps to use step sizes with absolute local discretization errors below `atol`. This tolerance must be a nonnegative numerical real value. The default tolerance is `AbsoluteError = 10^(-10*DIGITS)`.

The internal control mechanism estimates the local discretization error of a Runge-Kutta step and adjusts the step size adaptively to keep this error below the specified tolerances `rtol` or `atol`, respectively. The code uses the criterion

For accepting a solution vector Y. Roughly speaking, the relative error is controlled when the solution Y is sufficiently large. For very small solution values Y, absolute discretization errors are kept below the threshold `atol`.

Specify `AbsoluteError` = 0 if only control of the relative discretization errors is desired.

The error control may be switched off by specifying a fixed ```Stepsize = h```.

The default setting of , ensures that the local discretization errors are of the same order of magnitude as numerical roundoff.

Usually there is no need to use these options to change this setting. However, occasionally the numerical evaluation of the Runge-Kutta steps may be ill-conditioned or step sizes are so small that the time parameter cannot be incremented by the step size within working precision. In such a case these options may be used to bound the local discretization errors and use a higher working precision given by `DIGITS`.

Only positive real numerical values are accepted.

 Note:   The global error of the result returned by `numeric::odesolve` is usually larger than the local errors bounded by `rtol`, `atol`, respectively. Although the result is displayed with `DIGITS` decimal places one should not expect that all of them are correct. The relative precision of the final result is `rtol` at best!

See Example 8.

`Stepsize`

Option, specified as `Stepsize = h`

Switches off the internal error control and uses a Runge-Kutta iteration with constant step size `h` which must be a positive numerical value.

By default, `numeric::odesolve` uses an adaptive step size control mechanism in the numerical iteration. The option ```Stepsize = h``` switches off this adaptive mechanism and uses the Runge-Kutta method specified (or the default method `DOPRI78`) with constant step size `h`.

A final step with smaller step size is used to match the end `t` of the integration interval `t_0..t` if is not an integer.

 Note:   With this option, there is no automatic error control! Depending on the problem and on the order of the method the result may be a poor numerical approximation of the exact solution.

There is usually no need to invoke this option. However, occasionally the built-in adaptive error control mechanism may fail when integrating close to a singularity. In such a case this option may be used to customize a control mechanism for the global error by using different step sizes and investigating the convergence of the corresponding results.

Cf. Example 9.

`MaxStepsize`

Option, specified as `MaxStepsize = hmax`

Restricts adaptive step sizes to values not larger than `h_max`; `h_max` must be a positive numerical value.

By default, `numeric::odesolve` uses an adaptive step size control mechanism in the numerical iteration. The option ```MaxStepsize = hmax``` restricts the adaptive step size to values no larger than `hmax`.

If a larger stepsize `h` is requested explicitly by `Stepsize = h`, the option ```MaxStepsize = hmax``` reduces `h` to `hmax`.

`Alldata`

Option, specified as `Alldata = n`

Makes `numeric::odesolve` return a list of numerical mesh points generated by the internal Runge-Kutta iteration. The integer `n` controls the size of the output list.

With this option, `numeric::odesolve` returns a list of numerical mesh points [[t0, Y0], [t1, Y1], …, [t, Y(t)]] generated by the internal Runge-Kutta iteration.

The integer `n` controls the size of the output list. For n = 1, all internal mesh points are returned. This case may also be invoked by entering the simplified option `Alldata`, which is equivalent to `Alldata = 1`. For n > 1, only each n-th mesh point is stored in the list. The list always contains the initial point [t0, Y0] and the final point [t, Y(t)]. For n ≤ 0, only the data [[t0, Y0], [t, Y(t)]] are returned.

The output list may be useful to inspect the internal numerical process. Also further graphical processing of the mesh data may be useful.

Cf. Example 10.

`Symbolic`

Makes `numeric::odesolve` return a vector of symbolic expressions representing a single symbolic step of the Runge-Kutta iteration.

The call ```numeric::odesolve(f, t0..t, Y0, < method >, Symbolic)``` returns a vector (list or array) of expressions representing a single step of the numerical scheme with step size t - t0. In this mode symbolic values for t0, t and the components of Y0 are accepted.

This option may be useful if the specified numerical method applied to a given differential equation is to be investigated symbolically.

Cf. Example 11.

## Return Values

The solution vector Y(t) is returned as a list or as a 1-dimensional array of floating-point values. The type of the result vector coincides with the type of the input vector `Y0`.

With the option `Alldata`, a list of mesh data is returned.

## Algorithms

All methods presently implemented are adaptive versions of Runge-Kutta type single step schemes.

The methods `RKF43`, `RKF34`, `RKF54a`, `RKF54b`, `RKF45a`, `RKF45b`, `RKF87`, `RKF78`, `DOPRI54`, `DOPRI45`, `DOPRI65`, `DOPRI56`, `DOPRI87`, `DOPRI78`, `CK54`, `CK45` are embedded pairs of Runge-Kutta-Fehlberg, Dormand-Prince and Cash-Karp type, respectively. Estimates of the local discretization error are obtained in the usual way by comparing the results of the two submethods of the pair. The names indicate the orders of the subprocesses. For instance, `RKF34` and `RKF43` denote the same embedded Runge-Kutta-Fehlberg pair with orders 3 and 4. In `RKF34` the result of the fourth order submethod is accepted, whereas `RKF43` advances using the result of the third order submethod. In both cases the discretization error of the lower order subprocess is estimated and controlled.

For the single methods `EULER1` (the first order Euler scheme), `RK4` (the classical fourth order Runge-Kutta scheme) and `BUTCHER6` (a Runge-Kutta scheme of order 6), the relative local error is controlled by comparing steps with different step sizes. The effective order of these methods is increased by one through local extrapolation.

Local extrapolation is also available for the submethods of the embedded pairs. For instance, the method `xRKF78` uses extrapolation of the 8-th order subprocess of `RKF78`, yielding a method of effective order 9. The 7-th order subprocess is ignored. The cheap error estimate based on the embedded pair is not used implying some loss of efficiency when comparing `RKF78` (order 8) and `xRKF78` (effective order 9).

The call `numeric::butcher(method)` returns the Butcher data of the methods used in `numeric::odesolve`. Here `method` is one of `EULER1`, `RKF43`, `RK4`, `RKF34`, `RKF54a`, `RKF54b`, `DOPRI54`, `CK54`, `RKF45a`, `RKF45b`, `DOPRI45`, `CK45`, `DOPRI65`, `DOPRI56`, `BUTCHER6`, `RKF87`, `DOPRI87`, `RKF78`, `DOPRI78`.

 Note:   Only local errors are controlled by the adaptive mechanism. No control of the global error is provided!
 Note:   The run time of the numerical integration with a method of order p grows like , when `DIGITS` is increased. Computations with high precision goals are very expensive! High order methods such as the default method `DOPRI78` should be used.

Presently, only single step methods of Runge-Kutta type are implemented. Stiff problems cannot be handled efficiently with explicit methods such as the default method `DOPRI78`. For stiff problems, one may use one of the implicit A-stable methods `GAUSS(s)`. See Example 7.

For problems of the special type with a matrix valued function f(t, Y), there is a specialized ("geometric") integration routine `numeric::odesolveGeometric`. Generally, `numeric::odesolve` is faster than the geometric integrator. However, `odesolveGeometric` preserves certain invariants of the system more faithfully.

## References

J.C. Butcher: "The Numerical Analysis of Ordinary Differential Equations", Wiley, Chichester (1987).

E. Hairer, S.P. Norsett and G. Wanner: "Solving Ordinary Differential Equations I", Springer, Berlin (1993).