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

#### Introduction

`numeric::odesolve2``(f, t0, Y0, ..)` returns a function representing the numerical 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::odesolve2(f, t0, Y0 <, method> <, RelativeError = tol> <, Stepsize = h> )`

#### Parameters

`f` | - | procedure representing the vector field of the dynamical system
`t0` | - | numerical real value for the initial time.
`Y0` | - | list or 1-dimensional array of numerical values representing the initial value.

#### Options

`method` | - | name of the numerical scheme, see `numeric::odesolve`.
`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.

#### Returns

a procedure.

#### Side Effects

The function returned by `numeric::odesolve2` is sensitive to the environment variable `DIGITS`, which determines the numerical working precision. It uses `option remember`.

#### Related Functions

`numeric::odesolve`

#### Details

- The function generated by

  `Y := numeric::odesolve2(f, t0, Y0 <, options>)`

  is essentially

  `Y := t -> numeric::odesolve(t0..t, f, Y0 <,options>).`

- Numerical integration is launched, when `Y` is called with a real numerical argument. The call `Y(t)` returns the solution vector in a format corresponding to the type of the initial condition `Y0` with which `Y` was defined: `Y(t)` either yields a list or a 1-dimensional array.

  If `t` is not a real numerical value, then `Y(t)` returns an unevaluated function call.

- See the help page of `numeric::odesolve` for details on the parameters and the options.
- The options Alldata=`n` and Symbolic accepted by `numeric::odesolve` have no effect: `numeric::odesolve2` ignores these options.

- The function `Y` remembers all values it has computed. When calling `Y(t)` it searches its remember table for the time `T` closest to `t` and integrates from `T` to `t` using the previously computed `Y(T)` as initial value. This reduces the costs of a call considerably, if `Y` has to be evaluated many times (e.g., when plotting). However, the result `Y(t)` depends on the history of the MuPAD session! This can lead to unexpected side effects. Cf. example 3. We recommend to call `Y` only with a monotonically increasing (or decreasing) sequence of values `t` starting from `t0`. Further, the function must be re-initialized whenever `DIGITS` is increased. Cf. example 4.

- After the command `setuserinfo``(Y,1)` information on the current integration interval is displayed by each call to `Y`.

#### Example 1

The numerical solution of the initial value problem y'=t*sin(y), y(0)=2 is represented by the following function Y=[y]:

    >> f := (t, Y) -> [t*sin(Y[1])]: t0 := 0: Y0 := [2]:

    >> Y := numeric::odesolve2(f, 0, [2])

                                 proc Y(t) ... end

It starts numerical integration upon call with a numerical argument:

    >> Y(-2), Y(0), Y(0.1), Y(PI + sqrt(2))

                [2.968232567], [2.0], [2.004541745], [3.141552691]

Calling `Y` with a symbolic argument yields an unevaluated call:

    >> Y(t), Y(t + 5), Y(t^2 - 4)

                                                2
                             Y(t), Y(t + 5), Y(t  - 4)

    >> eval(subs(%, t = PI))

                    [3.13235701], [3.141592654], [3.141592611]

The numerical solution can be plotted. Note that `Y(t)` returns a list, so we plot the list element `Y(t)[1]`:

    >> plotfunc2d(Y(t)[1], t = -5..5)

    >> delete f, t0, Y0, Y:

#### Example 2

We consider 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=[Y1,Y2]=[y,y']:

    Y1'=Y2, Y2'=Y1^2, Y1(0)=0, Y2(0)=1.

    >> f := (t, Y) -> [Y[2], Y[1]^2]: t0 := 0: Y0 := [0, 1]:

    >> Y := numeric::odesolve2(f, t0, Y0):

    >> Y(1), Y(PI)

               [1.087473533, 1.362851121], [1274.867468, 37166.5226]

    >> delete f, t0, Y0, Y:

#### Example 3

We demonstrate a pitfall with the remember mechanism built into the functions returned by `numeric::odesolve2`. Consider the initial value problem y'=t*sin(y), y(0)=2:

    >> DIGITS := 5:
       Y := numeric::odesolve2((t, Y) -> [t*sin(Y[1])], 0, [2]):

The following result is stored in the remember table of `Y`:

    >> setuserinfo(Y, 1): Y(8.0)

          Info: integrating from t0=0 to t=8.0 using Y(t0)=[2]

                                     [3.1416]

The following value `Y(4.1)` is obtained using the previously computed `Y(8.0)`, integrating backward in time from t=8.0 to t=4.1:

    >> Y(4.1)

          Info: integrating from t0=8.0 to t=4.1 using Y(t0)=[3.1416]

                                     [4.1634]

We erase the remember table (the fifth operand) of `Y`:

    >> Y := subsop(Y, 5 = NIL):

The direct integration from the initial time t0=0 to t=4.1 yields a more accurate result than the previous computation:

    >> Y(4.1)

          Info: integrating from t0=0 to t=4.1 using Y(t0)=[2]

                                     [3.1413]

The reason for this phenomenon becomes obvious, when the solution of the ode is computed for various initial values Y(0)=2,3,4:

    >> A := numeric::odesolve2((t, Y) -> [t*sin(Y[1])], 0, [2]):
       B := numeric::odesolve2((t, Y) -> [t*sin(Y[1])], 0, [3]):
       C := numeric::odesolve2((t, Y) -> [t*sin(Y[1])], 0, [4]):

    >> plotfunc2d(A(t)[1], B(t)[1], C(t)[1], t = 0..8, Grid = 25)

In fact, all solutions with initial values Y(0) in the interval [0,2*PI] approach the same attracting point Y(infinity)=PI. While numerical integration forward in time is a very stable process, it is virtually impossible to recover the correct solution curve integrating backward in time starting at a point close to the attractor.

    >> setuserinfo(Y, 0): delete DIGITS, Y, A, B, C:

#### Example 4

We consider 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]]:
       Y := numeric::odesolve2(f, 0, [1, I]):

    >> DIGITS := 5: Y(1)

                      [3.5465 + 1.3683 I, 1.3683 + 0.80989 I]

Increasing `DIGITS` does not lead to a more accurate result because of the remember mechanism:

    >> DIGITS := 15: Y(1)

          [3.54647799893142 + 1.36829578207979 I,

             1.36829578207979 + 0.80988643477182 I]

This is the previous value computed with 5 digits, printed with 15 digits. Indeed, only 5 digits are correct. For getting a result that is accurate to full precision one has to erase the remember table via `Y:=subsop(Y,5=NIL)`. Alternatively, one may create a new numerical solution with a fresh (empty) remember table:

    >> Y := numeric::odesolve2(f, 0, [1, I]): Y(1)

          [3.54648242861716 + 1.36829887200859 I,

             1.36829887200859 + 0.80988468459998 I]

    >> delete f, Y, DIGITS:

#### Changes

- new function

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