Previous Page Next Page Contents

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

Option: method

Option: RelativeError = tol

Option: Stepsize = h

Option: Alldata = n

Option: Symbolic

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

Changes




Do you have questions or comments?


Copyright © SciFace Software GmbH & Co. KG 2000