## `numeric::odesolve` -- numerical solution of an ordinary differential equation

#### Introduction

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

           dY/dt = f(t,Y), Y(t0) = Y0, t0 and t real, Y and Y0 complex vectors.

#### Call(s)

`numeric::odesolve(t0..t, f, Y0 <, method> <, RelativeError = tol> <, Stepsize = h> <, Alldata = n> <, Symbolic> )`

#### Parameters

`t0` | - | numerical real value for the initial time
`t` | - | numerical real value (the ``time'')
`f` | - | procedure representing the vector field
`Y0` | - | list or 1-dimensional array of numerical values representing the initial condition Y0

#### Options

`method` | - | name of the numerical scheme
`RelativeError = tol` | - | forces internal numerical Runge-Kutta steps to use stepsizes with local discretization errors below `tol`. This tolerance must be a numerical real value >= 10^(-DIGITS).
`Stepsize = h` | - | switches off the internal error control and uses a Runge-Kutta iteration with constant stepsize `h`. `h` must be a numerical real value.
`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.
`Symbolic` | - | makes `numeric::odesolve` return a vector of symbolic expressions representing a single symbolic step of the Runge-Kutta iteration.

#### Returns

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.

#### Side Effects

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

#### Related Functions

`numeric::butcher`, `numeric::odesolve2`, `plot::ode`

#### Details

- `numeric::odesolve` is a general purpose solver able to deal with initial value problems of various kinds of (non-stiff) ordinary differential equations. Equations d^p/dt^p y(t)=f(t,y,y',..) of order p can be solved by `numeric::odesolve` after reformulation to dynamical system form. This can alway be achieved by writing the equation as a first order system

              d/dt [Y[1],..,Y[p-1],Y[p]] = [Y[2],..,Y[p],f(t,Y[1],..,Y[p])]

  for the vector [Y[1],Y[2],..]=[y,y',..]. Cf. example 4.
- 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 `exp(PI)`, `sqrt(2)` etc. are accepted.
- The vector field `f` defining the dynamical system Y'=f(t,Y) 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 `PI`, `exp(2)` etc. However, all values must be convertible to real or complex floating point numbers by `float`.
- Also autonomous systems, where f(t,Y) does not depend on t, must be represented by a procedure with 2 arguments `t` and `Y`.
- Also for scalar equations `Y` has to be represented by a list or an array with one element. For instance, the input data for the scalar initial value problem diff(y(t),t)=t*sin(y(t)), y(0)=1 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 stepsize keeps local relative discretization errors below `tol`=10^(-DIGITS), unless a different tolerance is specified via the option RelativeError=`tol`. The error control may be switched off by specifying a fixed Stepsize= `h`.

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

- With `Y := t -> numeric::odesolve(t0..t, f, Y0 <,options>)` 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 <, options>)`

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

- For systems of the special form diff(Y(t),t)=f(t,Y)*Y 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`.

#### Option: method

- Presently 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), | BUTCHER6 (order 6), | RKF87 (order 7),
  xRKF87 (order 7), | RKF78 (order 8), | xRKF78 (order 8).

  The Background below provides some information on these methods.

- There is usually no need to change the default method RKF78 by this option. This method is an embedded Runge-Kutta-Fehlberg pair of orders 7 and 8.
- Cf. example 6.

#### Option: RelativeError = tol

- The internal control mechanism estimates the local relative discretization error of a Runge-Kutta step and adjusts the stepsize adaptively to keep this error below `tol`.
- The default setting of `tol` = 10^(-DIGITS) ensures that the local discretization errors are of the same order of magnitude as numerical roundoff.

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

- Only real numerical values `tol` >= 10^(-DIGITS) are accepted.

- The global error of the result returned by `numeric::odesolve` is usually larger than the local errors bounded by `tol`. 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 `tol` at best!

- Cf. example 7.

#### Option: Stepsize = h

- By default, `numeric::odesolve` uses an adaptive stepsize control mechanism in the numerical iteration. The option Stepsize=`h` switches off this adaptive mechanism and uses the Runge-Kutta method specified by `method` (or the default method RKF78) with constant stepsize `h`.
- A final step with smaller stepsize is used to match the end `t` of the integration interval `t0..t`, if `(t-t0)/h` is not integer.

- 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 stepsizes and investigating the convergence of the corresponding results.
- Cf. example 8.

#### Option: Alldata = n

- 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 9.

#### Option: Symbolic

- The call `numeric::odesolve``(t0..t, f, Y0, <method>, Symbolic)` returns a vector (list or array) of expressions representing a single step of the numerical scheme with stepsize `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 10.

#### Example 1

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

    >> f := proc(t, Y) begin [t*sin(Y[1])] end_proc:

    >> numeric::odesolve(0..10, f, [2])

                                   [3.141592654]

    >> delete f:

#### Example 2

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

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

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

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

               [0.006737946999], [1.0], [2.718281828], [23.14069263]

    >> delete Y:

#### Example 3

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

     x'=x+y, y'=x-y, x(0)=1, y(0)=I.

    >> f := (t, Y) -> [Y[1] + Y[2], Y[1] - Y[2]]: Y0 := [1, I]:

    >> numeric::odesolve(0..PI, f, Y0)

            [72.57057162 + 30.05484302 I, 30.05484302 + 12.46088558 I]

The solution of a linear dynamical system Y'=A*Y with a constant matrix A is Y(t)=exp(t*A)*Y(0). 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)

            [72.57057163 + 30.05484303 I, 30.05484303 + 12.46088558 I]

    >> delete f, Y0, t, tA:

#### Example 4

We compute the numerical solution y(1) of the differential equation y''=y^2 with initial conditions y(0)=0, y'(0)=1 The second order equation is converted to a first order system for Y=[y,y']=[y,z]:

    y'=z, z'=y^2,y(0)=0, z(0)=1.

    >> f := proc(t, Y) begin [Y[2], Y[1]^2] end_proc: Y0 := [0, 1]:

    >> numeric::odesolve(0..1, f, Y0)

                            [1.087473533, 1.362851121]

    >> 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'=sin(t)-y, y(0)=1, the sample points are `t[0]=0.0`, `t[1]=0.1`, ..., `t[100]=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(t[i-1]..t[i], f, Y[i-1])
       end_for:

The following mesh data are produced:

    >> [t[i], Y[i]] $ i = 0..100

          [[0, [1]], [1/10, [0.9097]], [1/5, [0.8374]], [3/10, [0.7813]],

             [2/5, [0.7397]], ... , [99/10, [0.2159]], [10, [0.1476]]]

These data can be displayed by a list plot:

    >> plotpoints := [point(t[i], op(Y[i])) $ i = 0..100]:

    >> plot2d([Mode = List, plotpoints])

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

    >> plot(plot::ode(
           [t[i] $ i = 0..100], f, Y[0],
           [(t, Y) -> [t, Y[1]], Style = Points]))

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

#### Example 6

We compute the numerical solution y(1) of y'=y, y(0)=1 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(0..1, f, [1], RK4)

                                 [2.718281828459]

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(0..1, f, [1], xRKF78)

                                 [2.718281828459]

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 initial value problem y'=-10^(10)*y*(1-cos(y)), 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^10*Y[1]*(1 - cos(Y[1]))]: Y0 := [1]:

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

    >> DIGITS := 4: numeric::odesolve(0..1, f, Y0)

                                    [2.931e-5]

Due to numerical roundoff 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(0..1, f, Y0, RelativeError = 10^(-4))

                           [0.0000099999997577602193132]

Since relative local discretization errors are of the magnitude 10^(-4), 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(0..1, f, Y0, RelativeError = 10^(-10))

                           [0.000010000000000493745906]

    >> delete f, Y0, DIGITS:

#### Example 8

We compute the numerical solution y(1) of y'=y, y(0)=1 with various methods and various constant stepsizes. We compare the result with the exact solution y(1)=exp(1)=2.718281828....

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

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

    >> Y := numeric::odesolve(0..1, f, Y0, EULER1, Stepsize = 0.1):

    >> Y, globalerror = float(exp(1)) - Y[1]

                     [2.59374246], globalerror = 0.1245393684

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

    >> Y := numeric::odesolve(0..1, f, Y0, EULER1, Stepsize = 0.01):

    >> Y, globalerror = float(exp(1)) - Y[1]

                     [2.70481383], globalerror = 0.01346799904

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

    >> Y := numeric::odesolve(0..1, f, Y0, RK4, Stepsize = 0.1):

    >> Y, globalerror = float(exp(1)) - Y[1]

                  [2.718279744], globalerror = 0.000002084323879

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

    >> Y := numeric::odesolve(0..1, f, Y0, RK4, Stepsize = 0.01):

    >> Y, globalerror = float(exp(1)) - Y[1]

                [2.718281828], globalerror = 0.0000000002246438649

    >> delete f, Y0, Y:

#### Example 9

We integrate y'=y^2, y(0)=1 over the interval t=0..0.99 with a working precision of 4 digits. The exact solution is y(t)=1/(1-t). 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(0..0.99, f, Y0, Alldata)

          [[0.0, [1.0]], [0.5668, [2.308]], [0.7784, [4.513]],

             [0.8842, [8.636]], [0.9371, [15.9]], [0.9636, [27.43]],

             [0.9768, [43.05]], [0.99, [99.99]]]

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

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

          [[0.0, [1.0]], [0.7784, [4.513]], [0.9371, [15.9]],

             [0.9768, [43.05]], [0.99, [99.99]]]

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

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

          [[0.0, [1.0]], [0.1, [1.111]], [0.2, [1.25]], [0.3, [1.429]],

             [0.4, [1.667]], [0.5, [2.0]], [0.6, [2.5]], [0.7, [3.333]],

             [0.8, [5.0]], [0.9, [10.0]], [0.99, [131.2]]]

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

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

          [[0.0, [1.0]], [0.1, [1.111]], [0.2, [1.25]], [0.3, [1.429]],

             [0.4, [1.667]], [0.5, [2.0]], [0.6, [2.5]], [0.7, [3.333]],

             [0.8, [5.0]], [0.9, [10.0]], [0.99, [100.0]]]

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 10

The second order equation y''+sin(y)=0 is written as the dynamical system y'=z, z'=-sin(y) for the vector Y=[y,z]. A single symbolic step

    [y(t0),z(t0)] -> [y(t0+h), z(t0+h)]

of the Euler method is computed:

    >> f := proc(t, Y) begin [Y[2], -sin(Y[1])] end_proc:

    >> numeric::odesolve(t0..t0+h, f, [y0, z0], EULER1, Symbolic)

                            [y0 + h z0, z0 - h sin(y0)]

    >> delete f:

#### Background

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

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

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

- The methods RKF43, RKF34, RKF54a, RKF54b, RKF45a, RKF45b, RKF87, RKF78, DOPRI54, DOPRI45, 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 stepsizes. 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 RKF78 (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, BUTCHER6, RKF87, RKF78.

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

- The run time of the numerical integration with a method of order p grows like O(10^(DIGITS/(p+1))), when `DIGITS` is increased. Computations with high precision goals are very expensive! High order methods such as the default method RKF78 should be used.

- Presently only explicit single step methods of Runge-Kutta type are implemented. Stiff problems cannot be handled efficiently with these methods.
- For problems of the special type Y'=f(t,Y)*Y 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 sytem more faithfully.

#### Changes

- The new option RelativeError=`tol` allows to set a precision goal independent of the working precision.

Do you have questions or comments?
Copyright © SciFace Software GmbH & Co. KG 2000