Chapter 10

Ordinary Differential Equations

This chapter describes funtions for solving ordinary differential equation (ODE) initial value problems. The PLT Scheme Science Collection provides a variety of low-level methods, such as Runge-Kutta and Bulirsch-Stoer routines, and higher-level components for adaptive step-size control. The components can be combined by the user to achieve the desired solution, with full access to any intermediate steps. The functions described in this chapter are defined in the ode-initval.ss file in the science collection and are made available using the form:

(require (planet "ode-initval.ss" ("williams" "science.plt")))

10.1 Defining the ODE System

The routines solve the general n-dimensional first-order system,

dy_i(t)/dt = f_i(t, y_i(t), ... , y_n(t))

for i = 1, ..., n. The stepping functions rely on the vector of derivatives f_i and the Jacobian matrix, J_ij = df_i(t, y(t))/dy_j. A system of equations is defined using the ode-system structure.

ode-system

Structure:

ode-system

This structure defines a general ODE system with arbitrary parameters.


: function - A function (lambda (t y dydt params) ...). This function should store the elements of f_i(t, y, params) in the vector dydt, for arguments (t, y) and parameters params.
: jacobian - A function (lambda (t y dfdy dfdt params) .... This function should store the elements df_i(t, y, params)/dt in the vector dfdt and the Jacobian matrix J_ij in the vector dfdy as a row-ordered matrix J(i,j) = dfdy(i * dimension + j) where dimension is the dimension of the system. Some of the simplier solver algorithms do not make use of the Jacobian matrix, so it is not always strictly necessary to provide it (the jacobian field of the structure can be replace by #f for those algorithms). However, it is useful to provide the Jacobian to allow the solver algorithms to be interchanged. The best algorithms make use of the Jacobian.
: dimension - This is the dimension of the system of equations.
: params - This is a list of the arbitrary parameters of the system.

10.2 Stepping Functions

The lowest level components are the stepping functions that advance a solution from time t to t + h for a fixed step size h and estimate the resulting local error.

make-ode-step

Function:

(make-ode-step step-type dim)

Contract:

(-> ode-step-type? natural? ode-step?)

This function returns a newly created instance of a stepping function of type step-type for a system of dim dimensions.

ode-step-reset

Function:

(ode-step-reset step)

Contract:

(-> ode-step? void)

This function resets the stepping function step. It should be used whenever the next use of step will not be a continuation of a previous step.

ode-step-name

Function:

(ode-step-name step)

Contract:

(-> ode-step? string?)

This function returns the name of the stepping function step as a string.

ode-step-order

Function:

(ode-step-order step)

Contract:

(-> ode-step? natural?)

This function returns the order of the stepping function step on the previous step. This order can vary if the stepping function itself is adaptive.

ode-step-apply

Function:

(ode-step-apply step t h y y-err dydt-in dydt-out dydt)

Contract:

(-> ode-step? real? real?
    (vector-of real?) (vector-of real?)
    (vector-of real?) (vector-of real?)
    ode-system? void)

This function applies the stepping function step to the system of equations defined by dydt, using the step size h to advance the system from time t and state y to time t + h. The new state of the system is stored in y on output, with an estimate of the absolute error in each component stored in y-err. If the argument dydt-in is not #f, it should be a vector containing the derivatives for the system at time t on input. This is optional as the derivatives will be computed internally if they are not provided, but allows the reuse of existing derivative information. On output the new derivatives of the system at time t + h will be stored in the vector dydt-out, if it is not #f.

The following stepping algorithms are available.

Step Type:	`ode-step-rk2`

Embedded Runge-Kutta (2,3) method.

Step Type:	`ode-step-rk4`

4^th order (classical) Runge-Kutta.

Step Type:	`ode-step-rkf45`

Embedded Runge-Kutta-Fehlberg (4,5) method. This method is a good general-purpose integrator.

10.3 Adaptive Step-Size Control

The control function examines the proposed change to the solution and its error estimate produced by a stepping function and attempts to determine the optimal step-size for a user-specified level of error.

ode-control-standard-new

Function:

(ode-control-standard-new eps-abs eps-rel a-y a-dydt)

Contract:

(-> real? real? real? real?)

The standard control object is a four parameter heuristic based on absolute and relative errors eps-abs and eps-rel, and scaling factors a-y and a-dydt for the system stste y(t) and derivatives y'(t), respectively.

The step size adjustment procedure for this method begins by computing the desired error level D_i for each component,

$D_i = eps_{abs} + eps_{rel} * (a_y |y_i| + a_{dydt} h |y'_i|)$

and comparing it with the observed error E_i = |yerr_i|. If the observed error E exceeds the desired error level D by more than 10% for any component, then the method reduces the step size by an appropriate factor,

$h_{new} = h_{old} * S * (E/D)^{-1/q)}$

where q is the consistency order of the method (e.g. q=4 for 4(5) embedded RK), and S is a safety factor of 0.9. The ratio E/D is taken to be the maximum of the ratios E_i/D_i.

If the observed error E is less than 50% of the desired level D for the maximum ratio E_i/D_I, then the algorithm takes the opportunity to increase the step size to bring the error in line with the desired level,

$h_{new} = h_{old} * S * (E/D)^{(-1/(q+1))}.$

This encompasses all the standard scaling methods. To avoid uncontrolled changes in the step size, the overall scaling factor is limited to the range 1/5 to 5.

ode-control-y-new

Function:

(ode-control-y-new eps-abs eps-rel)

Contract:

(-> real? real?)

This function creates a new control object that will keep the local error within an absolute error of eps-abs and relative error eps-rel with respect to the solution y_i(t). This is equivalent to the standard control object with a-y = 1 and a-dydt = 0.

ode-control-yp-new

Function:

(ode-control-yp-new eps-abs eps-rel)

Contract:

(-> real? real?)

This function creates a new control object that will keep the local error within an absolute error of eps-abs and relative error eps-rel with respect to the derivatives of the solution y'_i(t). This is equivalent to the standard control object with a-y = 0 and a-dydt = 1.

ode-control-scaled-new

Function:

(ode-control-scaled-new eps-abs eps-rel a-y a-dydt
                        scale-abs dim)

Contract:

(-> real? real? real? real? (vector-of real?) natural?)

This function creates a new control object that uses the same algorithm as ode-control- standard-new, but with an absolute error that is scaled for each component by the array scale-abs. The formula for D_i for this control object is,

$D_i = eps_{abs} * s_i + eps_{rel} * (a_y |y_i| + a_{dydt} h |y'_i|)$

where s_i is the i^th component of the array scale-abs. The same error control heuristic is used by the Matlab ODE suite.

ode-control

Function:

(make-ode-control control-type)

Contract:

(-> ode-control-type? ode-control?)

This function returns a new instance of a control function of type control-type. This function is only needed for defining new types of control functions. For most purposes, the standard control functions described above should be sufficient.

ode-control-init

Function:

(ode-control-init control eps-abs eps-rel a-y a-dydt)

Contract:

(-> ode-control? real? real? real? real? any)

This function initializes the control function control with the parameters eps-abs (absolute error), eps-rel (relative error), a-y (scaling factor for y), and a-dydt scaling factor for derivatives.

ode-control-h-adjust

Function:

(ode-control-h-adjust control step y y-err dydt h)

Contract:

(-> ode-control? ode-step
    (vectorof real?) (vectorof real?)
     (vectorof real?) box? any)

This function adjusts the step size h using the control function control and the current values of y, y-err, and dydt. The stepping function step is also needed to determine the order of the method. If the error in the y valyes y-err is found to be too large, then the step-size h is reduced and the function returns - 1. If the error is sufficiently small, then h may be increased and 1 is returned. The function returns 0 if the step-size is unchanged. The goal of the function is to estimate the largest step-size that satisfies the user-specified accuracy requirement for the current point.

ode-control-name

Function:

(ode-control-name control)

Contract:

(-> ode-control? string?)

This function returns a pointer to the name of the control function. For eaxmple,

(printf ("control method is '~a'~n"
         (ode-control-name control)))

would print something like control mathod is 'standard'.

10.4 Evolution

The highest-level of the system is the evolution function that combines the results of a stepping function and control function to reliably advance the solution forward over an interval (t₀, t₁). If the control function signals that the step size should be descreased, the evolution function backs out of the current step and tries the proposed smaller step size. This process is continued until an acceptable step size is found.

make-ode-evolve

Function:

(make-ode-evolve dim)

Contract:

(-> positive? ode-evolve?)

This function returns a new instance of an evolution function for a system of dim dimensions.

ode-evolve-apply

Function:

(ode-evolve-apply evolve control step system
                  t t1 h y)

Contract:

(-> ode-evolve? ode-control? ode-step? ode-system?
    box? real? box? (vectorof real?) any)

This function evolves the system from time t and position y using the stepping function step. The new time and position are stored in t and y on output. The initial step size is taken as h, but this may be modified using the control function control to achieve the appropriate bound, if necessary. The routine may make several calls to step in order to determine the optimum step size. If the step size has been changed, the value of h will be modified on output. The maximum time t1 is guaranteed not to be exceeded by the time step. On the final time step, the value of t will be set to t1 exactly.

ode-evolve-reset

Function:

(ode-evolve-reset evolve)

Contract:

(-> ode-evolve? any)

This function resets the evolution function evolve. It should be used used whenever the next use of evolve will not be a continuation of a previous step.

10.5 Examples

Example: The following programs solve the second-order nonlinear Van der Pol oscillator equation,

This can be converted into a first order system suitable for use with the routines described in this chapter by introducing a separate variable for the velocity, y = x'(t),

$y' = -x + \mu y (1-x^2).$

The following example integrates the above system of equations from t = 0.0 to 100.0 in increments of 0.01 using a 4^th order Runge-Kutta stepping function.

(require (planet "ode-initval.ss" ("williams" "science.plt")))
(require (lib "plot.ss" "plot"))

(define (func t y f params)
  (let ((mu (car params))
        (y0 (vector-ref y 0))
        (y1 (vector-ref y 1)))
    (vector-set! f 0 y1)
    (vector-set! f 1 (- (- y0) (* mu y1 (- (* y0 y0) 1.0))))))

(define (main)
  (let* ((type rk4-ode-type)
         (step (make-ode-step type 2))
         (mu 10.0)
         (system (make-ode-system func #f 2 (list mu)))
         (t 0.0)
         (t1 100.0)
         (h 1.0e-2)
         (y #(1.0 0.0))
         (y-err (make-vector 2))
         (dydt-in (make-vector 2))
         (dydt-out (make-vector 2))
         (y0-values '())
         (y1-values '()))
    (ode-system-function-eval system t y dydt-in)
    (let loop ()
      (if (< t t1)
          (begin
            (ode-step-apply step t h
                               y y-err
                               dydt-in
                               dydt-out
                               system)
            ;(printf "~a ~a ~a~n" t (vector-ref y 0) (vector-ref y 1))
            (set! y0-values
                 (cons (vector t (vector-ref y 0)) y0-values))
            (set! y1-values
                 (cons (vector t (vector-ref y 1)) y1-values))
            (vector-set! dydt-in 0 (vector-ref dydt-out 0))
            (vector-set! dydt-in 1 (vector-ref dydt-out 1))
            (set! t (+ t h))
            (loop))))
    (printf "~a~n" (plot (points (reverse y0-values))
                         (x-min 0.0)
                         (x-max 100.0)
                         (y-min -2.0)
                         (y-max 2.0)))
    (printf "~a~n" (plot (points (reverse y1-values))
                         (x-min 0.0)
                         (x-max 100.0)))))

Figures 67 and 68 show the resulting output plots of y0 and y1.

Figure 67: ODE Example 1 Plot of y0

Figure 68: ODE Example 1 Plot of y1

Example: The following example evolves the above system of equations from t = 0.0 to 100.0 maintaining an error in the y value of 1.0e-6 using a 4^th order Runge-Kutta stepping function.

(require (planet "ode-initval.ss" ("williams" "science.plt")))
(require (lib "plot.ss" "plot"))

(define (func t y f params)
  (let ((mu (car params))
        (y0 (vector-ref y 0))
        (y1 (vector-ref y 1)))
    (vector-set! f 0 y1)
    (vector-set! f 1 (- (- y0) (* mu y1 (- (* y0 y0) 1.0))))))

(define (main)
  (let* ((type rk4-ode-type)
         (step (make-ode-step type 2))
         (control (control-y-new 1.0e-6 0.0))
         (evolve (make-ode-evolve 2))
         (mu 10.0)
         (system (make-ode-system func #f 2 (list mu)))
         (t (box 0.0))
         (t1 100.0)
         (h (box 1.0e-6))
         (y #(1.0 0.0))
         (y0-values '())
         (y1-values '()))
    (let loop ()
      (if (< (unbox t) t1)
          (begin
            (ode-evolve-apply
             evolve control step system
             t t1 h y)
            ;(printf "~a ~a ~a~n"
            ;       (unbox t) (vector-ref y 0) (vector-ref y 1))
            (set! y0-values
                 (cons (vector (unbox t) (vector-ref y 0)) y0-values))
            (set! y1-values
                 (cons (vector (unbox t) (vector-ref y 1)) y1-values))
            (loop))))
    (printf "Number of iterations   = ~a~n"
            (ode-evolve-count evolve))
    (printf "Number of failed steps = ~a~n"
            (ode-evolve-failed-steps evolve))
    (printf "~a~n" (plot (points (reverse y0-values))
                         (x-min 0.0)
                         (x-max 100.0)
                         (y-min -2.0)
                         (y-max 2.0)))
    (printf "~a~n" (plot (points (reverse y1-values))
                         (x-min 0.0)
                         (x-max 100.0)))))

When run, it prints the following:

Number of iterations   = 84575
Number of failed steps = 352

Figures 69 and 70 show the resulting output plots of y0 and y1.

Figure 69: ODE Example 2 Plot of y0

Figure 70: ODE Example 2 Plot of y1

Chapter 10 Ordinary Differential Equations

10.1 Defining the ODE System

10.2 Stepping Functions

10.3 Adaptive Step-Size Control

10.4 Evolution

10.5 Examples

Chapter 10

Ordinary Differential Equations