scheme-example.texinfo (3523B)
1 @node Example 2 @chapter Example 3 4 The procedure @code{integrate-system} integrates the system 5 6 @displaymath 7 y_k^\prime = f_k(y_1, y_2, \ldots, y_n), \; k = 1, \ldots, n 8 @end displaymath 9 10 of differential equations with the method of Runge-Kutta. 11 12 The parameter @code{system-derivative} is a function that takes a 13 system state (a vector of values for the state variables 14 @math{y_1, \ldots, y_n})and produces a system derivative (the values 15 @math{y_1^\prime, \ldots, y_n^\prime}). The parameter 16 @code{initial-state} provides an initial system state, and @code{h} 17 is an initial guess for the length of the integration step. 18 19 The value returned by @code{integrate-system} is an infinite stream of 20 system states. 21 22 @lisp 23 (define (integrate-system system-derivative 24 initial-state 25 h) 26 (let ((next (runge-kutta-4 system-derivative h))) 27 (letrec ((states 28 (cons initial-state 29 (delay (map-streams next 30 states))))) 31 states))) 32 @end lisp 33 34 The procedure @code{runge-kutta-4} takes a function, @code{f}, that 35 produces a system derivative from a system state. It produces a 36 function that takes a system state and produces a new system state. 37 38 @lisp 39 (define (runge-kutta-4 f h) 40 (let ((*h (scale-vector h)) 41 (*2 (scale-vector 2)) 42 (*1/2 (scale-vector (/ 1 2))) 43 (*1/6 (scale-vector (/ 1 6)))) 44 (lambda (y) 45 ;; y is a system state 46 (let* ((k0 (*h (f y))) 47 (k1 (*h (f (add-vectors y (*1/2 k0))))) 48 (k2 (*h (f (add-vectors y (*1/2 k1))))) 49 (k3 (*h (f (add-vectors y k2))))) 50 (add-vectors y 51 (*1/6 (add-vectors k0 52 (*2 k1) 53 (*2 k2) 54 k3))))))) 55 56 (define (elementwise f) 57 (lambda vectors 58 (generate-vector 59 (vector-length (car vectors)) 60 (lambda (i) 61 (apply f 62 (map (lambda (v) (vector-ref v i)) 63 vectors)))))) 64 65 (define (generate-vector size proc) 66 (let ((ans (make-vector size))) 67 (letrec ((loop 68 (lambda (i) 69 (cond ((= i size) ans) 70 (else 71 (vector-set! ans i (proc i)) 72 (loop (+ i 1))))))) 73 (loop 0)))) 74 75 (define add-vectors (elementwise +)) 76 77 (define (scale-vector s) 78 (elementwise (lambda (x) (* x s)))) 79 @end lisp 80 81 The @code{map-streams} procedure is analogous to @code{map}: it applies 82 its first argument (a procedure) to all the elements of its second 83 argument (a stream). 84 85 @lisp 86 (define (map-streams f s) 87 (cons (f (head s)) 88 (delay (map-streams f (tail s))))) 89 @end lisp 90 91 Infinite streams are implemented as pairs whose car holds the first 92 element of the stream and whose cdr holds a promise to deliver the rest 93 of the stream. 94 95 @lisp 96 (define head car) 97 (define (tail stream) 98 (force (cdr stream))) 99 @end lisp 100 101 The following illustrates the use of @code{integrate-system} in 102 integrating the system 103 104 @displaymath 105 C {dv_C \over dt} = -i_L - {v_C \over R} 106 @end displaymath 107 108 @displaymath 109 L {di_L \over dt} = v_C 110 @end displaymath 111 112 which models a damped oscillator. 113 114 @lisp 115 (define (damped-oscillator R L C) 116 (lambda (state) 117 (let ((Vc (vector-ref state 0)) 118 (Il (vector-ref state 1))) 119 (vector (- 0 (+ (/ Vc (* R C)) (/ Il C))) 120 (/ Vc L))))) 121 122 (define the-states 123 (integrate-system 124 (damped-oscillator 10000 1000 .001) 125 '#(1 0) 126 .01)) 127 @end lisp