例

integrate-system : 次の微文方程式系

  y'_k = f_k(y_1, y_2, ..., y_n), k = 1, ..., n

を、ルンゲ・クッタ(Runge-Kutta)法を使用して積分する。

パラメータsystem-derivativeは、系の状態(状態変数y_1, ...,y_nについての値のベクタ)を引数に取って、系の導関数(値y'_1, ...,y'_n)を生成する。パラメータinitial-stateは、系の初期状態を与える。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
                       (delay (map-streams next
                                        states))))) 
        states))))

Runge-Kutta-4は、系の状態から系の導関数を生成する関数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は系の状態
        (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に似ており、第1引数(手続き)を第2引数(ストリーム)のすべての要素に適用する。

(define map-streams
  (lambda (f s)
    (cons (f (head s))
          (delay (map-streams f (tail s))))))

無限ストリームはペアとして実装され、そのcarにはストリームの最初の要素が、cdrにはストリームの残りを生成するプロミスが保持される。

(define head car) 
(define tail
  (lambda (stream) (force (cdr stream))))

integrate-systemの利用法を下記に示す。これは減衰発信器をモデル化した次の系の積分を行なうものである。

  C * (dv_C / dt) = -i_L - (v_C / R)
  L * (di_L / dt) = v_C

(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)