本節では (runge-kutta) ライブラリから成る例を説明する。これは微分方程式の系
y'_k = f_k(y_1, y_2, ..., y_n), k = 1, ..., n
を Runge-Kutta 方で積分する積分系手続きを提供するものである。
(runge-kutta) ライブラリは (rnrs base (6)) を使うため、その概形は以下のようになる。
#!r6rs (library (runge-kutta) (export integrate-system head tail) (import (rnrs base)) <library body>)
以下で説明する手続きの定義は <library body> の部分に入る。
仮引き数 system-derivative は系の状態(状態変数 y_1, ..., y_n のベクタ)を受け取り、系の微分係数(y '_1, ..., y '_n)を返す。仮引き数は初期状態を与え、 h は積分回数の初期的な推測である。
integrate-system の戻り値は系の状態の無限ストリームである。
(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))))
runge-kutta-4 手続きは関数 f を取り、 f は系の状態から系の微分係数を算出する。 runge-kutta-4 手続きは系の状態を受け取り系の新しい状態を返す関数を返す。
(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 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)))))
map-streams 手続きは map に類似のものである。最初の引き数(手続き)を 2 番目の手続き(ストリーム)の要素すべてに適用する。
(define map-streams (lambda (f s) (cons (f (head s)) (lambda () (map-streams f (tail s))))))
無限ストリームは car がストリームの最初の要素であり、 cdr がストリームの残りの要素を返す手続きを保持する対として実装されている。
(define head car) (define tail (lambda (stream) ((cdr stream))))
以下のプログラムは次の系を積分する integrating-system の使用法を説明している。
C (dv_C / dt) = -i_L (-v_C / R) L (di_L / dt) = v_C
これは減衰振動をモデル化したものである。
#!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))
このプログラムは次のようなものを表示する。
#(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)