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 lowlevel methods, such as RungeKutta and BulirschStoer routines, and higherlevel components for adaptive stepsize 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 odeinitval.ss file in the science collection and are made available using the form:
(require (planet "odeinitval.ss" ("williams" "science.plt")))
10.1 Defining the ODE System
The routines solve the general ndimensional firstorder system,
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 odesystem structure.
Structure:
odesystem 

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.
Function:
(makeodestep steptype dim) 
Contract: (> odesteptype? natural? odestep?) 

Function:
(odestepreset step) 
Contract: (> odestep? void) 

Function:
(odestepname step) 
Contract: (> odestep? string?) 

Function:
(odesteporder step) 
Contract: (> odestep? natural?) 

Function:
(odestepapply step t h y yerr dydtin dydtout dydt) 
Contract: (> odestep? real? real? (vectorof real?) (vectorof real?) (vectorof real?) (vectorof real?) odesystem? void) 

The following stepping algorithms are available.
Embedded RungeKutta (2,3) method.
4^{th} order (classical) RungeKutta.
Embedded RungeKuttaFehlberg (4,5) method. This method is a good generalpurpose integrator.
10.3 Adaptive StepSize 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 stepsize for a userspecified level of error.
Function:
(odecontrolstandardnew epsabs epsrel ay adydt) 
Contract: (> real? real? real? real?) 

The step size adjustment procedure for this method begins by computing the desired error level D_{i} for each component,
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,
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,
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.
Function:
(odecontrolynew epsabs epsrel) 
Contract: (> real? real?) 

Function:
(odecontrolypnew epsabs epsrel) 
Contract: (> real? real?) 

Function:
(odecontrolscalednew epsabs epsrel ay adydt scaleabs dim) 
Contract: (> real? real? real? real? (vectorof real?) natural?) 

where s_{i} is the i^{th} component of the array scaleabs. The same error control heuristic is used by the Matlab ODE suite.
Function:
(makeodecontrol controltype) 
Contract: (> odecontroltype? odecontrol?) 

Function:
(odecontrolinit control epsabs epsrel ay adydt) 
Contract: (> odecontrol? real? real? real? real? any) 

Function:
(odecontrolhadjust control step y yerr dydt h) 
Contract: (> odecontrol? odestep (vectorof real?) (vectorof real?) (vectorof real?) box? any) 

Function:
(odecontrolname control) 
Contract: (> odecontrol? string?) 

(printf ("control method is '~a'~n" (odecontrolname control)))
would print something like control mathod is 'standard'.
10.4 Evolution
The highestlevel 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_{0}, t_{1}). 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.
Function:
(makeodeevolve dim) 
Contract: (> positive? odeevolve?) 

Function:
(odeevolveapply evolve control step system t t1 h y) 
Contract: (> odeevolve? odecontrol? odestep? odesystem? box? real? box? (vectorof real?) any) 

Function:
(odeevolvereset evolve) 
Contract: (> odeevolve? any) 

10.5 Examples
Example: The following programs solve the secondorder 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),
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 RungeKutta stepping function.
(require (planet "odeinitval.ss" ("williams" "science.plt"))) (require (lib "plot.ss" "plot")) (define (func t y f params) (let ((mu (car params)) (y0 (vectorref y 0)) (y1 (vectorref y 1))) (vectorset! f 0 y1) (vectorset! f 1 ( ( y0) (* mu y1 ( (* y0 y0) 1.0)))))) (define (main) (let* ((type rk4odetype) (step (makeodestep type 2)) (mu 10.0) (system (makeodesystem func #f 2 (list mu))) (t 0.0) (t1 100.0) (h 1.0e2) (y #(1.0 0.0)) (yerr (makevector 2)) (dydtin (makevector 2)) (dydtout (makevector 2)) (y0values '()) (y1values '())) (odesystemfunctioneval system t y dydtin) (let loop () (if (< t t1) (begin (odestepapply step t h y yerr dydtin dydtout system) ;(printf "~a ~a ~a~n" t (vectorref y 0) (vectorref y 1)) (set! y0values (cons (vector t (vectorref y 0)) y0values)) (set! y1values (cons (vector t (vectorref y 1)) y1values)) (vectorset! dydtin 0 (vectorref dydtout 0)) (vectorset! dydtin 1 (vectorref dydtout 1)) (set! t (+ t h)) (loop)))) (printf "~a~n" (plot (points (reverse y0values)) (xmin 0.0) (xmax 100.0) (ymin 2.0) (ymax 2.0))) (printf "~a~n" (plot (points (reverse y1values)) (xmin 0.0) (xmax 100.0)))))
Figures 67 and 68 show the resulting output plots of y0 and 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.0e6 using a 4^{th} order RungeKutta stepping function.
(require (planet "odeinitval.ss" ("williams" "science.plt"))) (require (lib "plot.ss" "plot")) (define (func t y f params) (let ((mu (car params)) (y0 (vectorref y 0)) (y1 (vectorref y 1))) (vectorset! f 0 y1) (vectorset! f 1 ( ( y0) (* mu y1 ( (* y0 y0) 1.0)))))) (define (main) (let* ((type rk4odetype) (step (makeodestep type 2)) (control (controlynew 1.0e6 0.0)) (evolve (makeodeevolve 2)) (mu 10.0) (system (makeodesystem func #f 2 (list mu))) (t (box 0.0)) (t1 100.0) (h (box 1.0e6)) (y #(1.0 0.0)) (y0values '()) (y1values '())) (let loop () (if (< (unbox t) t1) (begin (odeevolveapply evolve control step system t t1 h y) ;(printf "~a ~a ~a~n" ; (unbox t) (vectorref y 0) (vectorref y 1)) (set! y0values (cons (vector (unbox t) (vectorref y 0)) y0values)) (set! y1values (cons (vector (unbox t) (vectorref y 1)) y1values)) (loop)))) (printf "Number of iterations = ~a~n" (odeevolvecount evolve)) (printf "Number of failed steps = ~a~n" (odeevolvefailedsteps evolve)) (printf "~a~n" (plot (points (reverse y0values)) (xmin 0.0) (xmax 100.0) (ymin 2.0) (ymax 2.0))) (printf "~a~n" (plot (points (reverse y1values)) (xmin 0.0) (xmax 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.

