diff options
| -rw-r--r-- | test/damped-harmonic.scm | 78 |
1 files changed, 78 insertions, 0 deletions
diff --git a/test/damped-harmonic.scm b/test/damped-harmonic.scm new file mode 100644 index 0000000..a47be34 --- /dev/null +++ b/test/damped-harmonic.scm @@ -0,0 +1,78 @@ + +; This is the "Example" from R7RS-little + +(define (integrate-system system-derivative + initial-state + h) + (let ((next (runge-kutta-4 system-derivative h))) + (letrec ((states + (cons initial-state + (delay (map-streams next + states))))) + states))) + +(define (runge-kutta-4 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 f) + (lambda vectors + (generate-vector + (vector-length (car vectors)) + (lambda (i) + (apply f + (map (lambda (v) (vector-ref v i)) + vectors)))))) + +(define (generate-vector size proc) + (let ((ans (make-vector size))) + (letrec ((loop + (lambda (i) + (cond ((= i size) ans) + (else + (vector-set! ans i (proc i)) + (loop (+ 1 i))))))) + (loop 0)))) + +(define add-vectors (elementwise +)) + +(define (scale-vector s) + (elementwise (lambda (x) (* x s)))) + +(define (map-streams f s) + (cons (f (head s)) + (delay (map-streams f (tail s))))) + +(define head car) +(define (tail stream) + (force (cdr stream))) + +; Ok, now the actual example + +(define (damped-oscillator 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)) + +(display the-states) |