example.tex

\chapter{Example }
\label{exampleappendix}

\nobreak
This section describes an example consisting of the
\library{runge-kutta} library, which provides an {\cf integrate-system}
procedure that integrates the system 
$$y_k^\prime = f_k(y_1, y_2, \ldots, y_n), \; k = 1, \ldots, n$$
of differential equations with the method of Runge-Kutta.

As the \library{runge-kutta} library makes use of the \rsixlibrary{base}
library, its skeleton is as follows:

\begin{scheme}
\#!r6rs
(library (runge-kutta)
  (export integrate-system
          head tail)
  (import (rnrs base))
  \hyper{library body})
\end{scheme}

The procedure definitions described below go in the place of \hyper{library body}.

The parameter {\tt system-derivative} is a function that takes a system
state (a vector of values for the state variables $y_1, \ldots, y_n$)
and produces a system derivative (the values $y_1^\prime, \ldots,
y_n^\prime$).  The parameter {\tt initial-state} provides an initial
system state, and {\tt h} is an initial guess for the length of the
integration step.

The value returned by {\cf integrate-system} is an infinite stream of
system states.

\begin{schemenoindent}
(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))))%
\end{schemenoindent}

The {\cf runge-kutta-4} procedure takes a function, {\tt f}, that produces a
system derivative from a system state.  The {\cf runge-kutta-4} procedure
produces a function that takes a system state and
produces a new system state.

\begin{schemenoindent}
(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 {\rm{}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)))))%
\end{schemenoindent}

The {\cf map-streams} procedure is analogous to {\cf map}: it applies its first
argument (a procedure) to all the elements of its second argument (a
stream).

\begin{schemenoindent}
(define map-streams
  (lambda (f s)
    (cons (f (head s))
          (lambda () (map-streams f (tail s))))))%
\end{schemenoindent}

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.

\begin{schemenoindent}
(define head car)
(define tail
  (lambda (stream) ((cdr stream))))%
\end{schemenoindent}

\bigskip
The following program illustrates the use of {\cf integrate-\hp{}system} in
integrating the system
$$ C {dv_C \over dt} = -i_L - {v_C \over R}$$\nobreak
$$ L {di_L \over dt} = v_C$$
which models a damped oscillator.

\begin{schemenoindent}
\#!r6rs
(import (rnrs base)
        (rnrs io simple)
        (runge-kutta))

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

(letrec ((loop (lambda (s)
                 (newline)
                 (write (head s))
                 (loop (tail s)))))
  (loop the-states))%
\end{schemenoindent}

This prints output like the following:

\begin{scheme}
\#(1 0)
\#(0.99895054 9.994835e-6)
\#(0.99780226 1.9978681e-5)
\#(0.9965554 2.9950552e-5)
\#(0.9952102 3.990946e-5)
\#(0.99376684 4.985443e-5)
\#(0.99222565 5.9784474e-5)
\#(0.9905868 6.969862e-5)
\#(0.9888506 7.9595884e-5)
\#(0.9870173 8.94753e-5)
\end{scheme}

%%% Local Variables: 
%%% mode: latex
%%% TeX-master: "r6rs"
%%% End: