This section describes an example consisting of the
(runge-kutta) library, which provides an integrate-system
procedure that integrates the system

yk⁄ = fk(y1, y2, ..., yn), ‌k = 1, ..., n

of differential equations with the method of Runge-Kutta.

As the (runge-kutta) library makes use of the (rnrs base (6))
library, its skeleton is as follows:

#!r6rs

(library (runge-kutta)

(export integrate-system

head tail)

(import (rnrs base))

<library body>)

The procedure definitions described below go in the place of <library body>.

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 y1⁄, ...,
yn⁄). 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

(lambda ()

(map-streams next states)))))

states))))

The runge-kutta-4 procedure takes a function, f, that produces a
system derivative from a system state. The runge-kutta-4 procedure
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)))))

The map-streams procedure 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))

(lambda () (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 procedure that delivers the rest
of the stream.

(define head car)

(define tail

(lambda (stream) ((cdr stream))))

The following program illustrates the use of integrate-system in
integrating the system