R5Rs
Scheme Revised(5) Report on the Algorithmic Language Scheme -- Example
Integrate-system integrates the system

y'k = fk(y1, y2, ..., yn), k = 1, ..., n
of differential equations with the method of Runge-Kutta.

The parameter system-derivative is a function that takes a system state (a vector of values for the state variables y1, ..., yn) and produces a system derivative (the values y'1, ..., y'n). The parameter initial-state provides an initial system state, and h is an initial guess for the length of the integration step.

The value returned by integrate-system is an infinite stream of system states.

(define integrate-system
  (lambda (system-derivative initial-state h)
    (let ((next (runge-kutta-4 system-derivative h)))
      (letrec ((states
                (cons initial-state
                      (delay (map-streams next
                                          states)))))
        states))))
Runge-Kutta-4 takes a function, f, that produces a system derivative from a system state. Runge-Kutta-4 produces a function that takes a system state and produces a new system state.


(define runge-kutta-4
  (lambda (f h)
    (let ((*h (scale-vector h))
          (*2 (scale-vector 2))
          (*1/2 (scale-vector (/ 1 2)))
          (*1/6 (scale-vector (/ 1 6))))
      (lambda (y)
        ;; y is a system state
        (let* ((k0 (*h (f y)))
               (k1 (*h (f (add-vectors y (*1/2 k0)))))
               (k2 (*h (f (add-vectors y (*1/2 k1)))))
               (k3 (*h (f (add-vectors y k2)))))
          (add-vectors y
            (*1/6 (add-vectors k0
                               (*2 k1)
                               (*2 k2)
                               k3))))))))

(define elementwise (lambda (f) (lambda vectors (generate-vector (vector-length (car vectors)) (lambda (i) (apply f (map (lambda (v) (vector-ref v i)) vectors)))))))

(define generate-vector (lambda (size proc) (let ((ans (make-vector size))) (letrec ((loop (lambda (i) (cond ((= i size) ans) (else (vector-set! ans i (proc i)) (loop (+ i 1))))))) (loop 0)))))

(define add-vectors (elementwise +))

(define scale-vector (lambda (s) (elementwise (lambda (x) (* x s)))))
Map-streams is analogous to map: it applies its first argument (a procedure) to all the elements of its second argument (a stream).

(define map-streams
  (lambda (f s)
    (cons (f (head s))
          (delay (map-streams f (tail s))))))
Infinite streams are implemented as pairs whose car holds the first element of the stream and whose cdr holds a promise to deliver the rest of the stream.

(define head car)
(define tail
  (lambda (stream) (force (cdr stream))))
The following illustrates the use of integrate-system in integrating the system

C dvC / dt = -iL - vC / R
L diL / dt = vC
which models a damped oscillator.

(define damped-oscillator
  (lambda (R L C)
    (lambda (state)
      (let ((Vc (vector-ref state 0))
            (Il (vector-ref state 1)))
        (vector (- 0 (+ (/ Vc (* R C)) (/ Il C)))
                (/ Vc L))))))

(define the-states (integrate-system (damped-oscillator 10000 1000 .001) '#(1 0) .01))

This Html page has been produced by Skribe.
Last update Sat Jul 13 09:52:06 2013.