R6RS:翻訳:R6RS:D Example

R6RS:翻訳:R6RS:D Example

付録 D 例

本節では (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)

Last modified : 2008/05/05 10:35:14 UTC