## Example

‘`Integrate-system`’ integrates the system

y_k^^ = f_k(y_1, y_2, ..., y_n), 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 y_1, ..., y_n) 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 dv_C / dt = -i_L - v_C / R
L di_L / dt = v_C

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))