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