Previous Page Next Page Contents

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

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




Do you have questions or comments?


Copyright © SciFace Software GmbH & Co. KG 2000