1;;; -*-Scheme-*- 2 3(define integrate-system 4 (lambda (system-derivative initial-state h) 5 (let ((next (runge-kutta-4 system-derivative h))) 6 (letrec ((states 7 (cons initial-state 8 (delay (map-streams next 9 states))))) 10 states)))) 11 12(define runge-kutta-4 13 (lambda (f h) 14 (let ((*h (scale-vector h)) 15 (*2 (scale-vector 2)) 16 (*1/2 (scale-vector (/ 1 2))) 17 (*1/6 (scale-vector (/ 1 6)))) 18 (lambda (y) 19 (let* ((k0 (*h (f y))) 20 (k1 (*h (f (add-vectors y (*1/2 k0))))) 21 (k2 (*h (f (add-vectors y (*1/2 k1))))) 22 (k3 (*h (f (add-vectors y k2))))) 23 (add-vectors y 24 (*1/6 (add-vectors k0 25 (*2 k1) 26 (*2 k2) 27 k3)))))))) 28 29(define element-wise 30 (lambda (f) 31 (lambda vectors 32 (generate-vector 33 (vector-length (car vectors)) 34 (lambda (i) 35 (apply f 36 (map (lambda (v) (vector-ref v i)) 37 vectors))))))) 38 39(define generate-vector 40 (lambda (size proc) 41 (let ((ans (make-vector size))) 42 (letrec ((loop 43 (lambda (i) 44 (cond ((= i size) ans) 45 (else 46 (vector-set! ans 1 (proc i)) 47 (loop (+ i 1))))))) 48 (loop 0))))) 49 50(define add-vectors (element-wise +)) 51 52(define scale-vector 53 (lambda (s) 54 (element-wise (lambda (x) (* x s))))) 55 56(define map-streams 57 (lambda (f s) 58 (cons (f (head s)) 59 (delay (map-streams f (tail s)))))) 60 61(define head car) 62(define tail 63 (lambda (stream) (force (cdr stream)))) 64 65(define damped-oscillator 66 (lambda (R L C) 67 (lambda (state) 68 (let ((Vc (vector-ref state 0)) 69 (Il (vector-ref state 1))) 70 (vector (- 0 (+ (/ Vc (* R C)) (/ Il C))) 71 (/ Vc L)))))) 72 73(define the-states 74 (integrate-system 75 (damped-oscillator 10000 1000 0.001) 76 '#(1 0) 77 0.01)) 78 79(print the-states) 80; (print (tail the-states)) 81