SICP 読書ノート#39 - 3.5.4 ストリームと遅延評価 (pp.205-209)

(2015/09/15追記) SICPテキストのストリームでは問題3.77のテスト実行で返ってこないため、ここではracket/stream版を使っています。ソースコードはGitHubに置いています。

前節で出てきた積分器integralを使ってフィードバックループを持つ信号処理システムを考える。大学では制御工学を専攻していたのでこういうの大好きです。

このフィードバックシステムは以下のように実装される。

(define (solve f y0 dt)
  (define y (integral (delay dy) y0 dt))
  (define dy (stream-map f y))
  y)

積分器integralの入力ストリームをdelayしているところがミソ。現在時刻をtとするとdy(t-1)を入力する必要があるためこれで良い。

積分器は入力ストリームをforceで評価する必要がある。

(define (integral delayed-integrand initial-value dt)
  (define int
    (cons-stream initial-value
                 (let ((integrand (force delayed-integrand)))
                   (add-streams (scale-stream integrand dt)
                                int))))
  int)

テスト。

racket@> (stream-ref (solve (lambda (y) y) 1 0.001) 1000)
=> 2.716923932235896

たぶんdtの値を小さくするほど \( e \) に近づくはず。

問題 3.77

テキストで定義されているintegers-starting-fromに似たintegral手続きを遅延リストに対応させよ。

(define (integral-ex delayed-integrand initial-value dt)
  (cons-stream initial-value
               (let ((integrand (force delayed-integrand)))
                 (if (stream-null? integrand)
                     the-empty-stream
                     (integral-ex (stream-cdr integrand)
                                  (+ initial-value
                                     (* (stream-car integrand) dt))
                                  dt)))))

(define (solve-ex f y0 dt)
  (define y (integral-ex (delay dy) y0 dt))
  (define dy (stream-map f y))
  y)

テスト。先程の結果と同様。

(define (solve-ex f y0 dt)
  (define y (integral-ex (delay dy) y0 dt))
  (define dy (stream-map f y))
  y)
;;=> 2.716923932235896

問題 3.78

以下の二次微分方程式

\( \frac{d^2 y}{dt^2} - a \frac{dy}{dt} - by = 0\)

の解を求めるためのフィードバックシステムで、yを求める手続きsolve-2ndを実装せよ。

テキストの図のフィードバックシステムの通りに実装すればよいので、

(define (solve-2nd dy0 y0 a b dt)
  (define y (integral (delay dy) y0 dt))
  (define dy (integral
              (delay (add-streams
                      (scale-stream dy a)
                      (scale-stream y b)))
              dy0 dt))
  y)

テスト。

racket@> (stream-ref (solve-2nd 1 1 0 1 0.001) 1000)
=> 2.716923932235896
racket@> (stream-ref (solve-2nd 1 1 2 -1 0.001) 1000)
=> 2.716923932235896

問題 3.79

問題3.78をさらに汎用化して

\( \frac{d^2 y}{dt^2} = f( \frac{dy}{dt} , y ) \)

のフィードバックシステムにおいて、yを求める手続きを実装する。

(define (solve-2nd-ex f y0 dy0 dt)
  (define y (integral (delay dy) y0 dt))
  (define dy (integral (delay ddy) dy0 dt))
  (define ddy (stream-map f dy y))
  y)

テスト。fを問題3.78と同じものを与えれば、当然結果も同じになる。

racket@> (stream-ref
          (solve-2nd-ex (lambda (dy y) y)
                        1 1 0.001) 1000)
=> 2.716923932235896
racket@> (stream-ref
          (solve-2nd-ex (lambda (dy y) (+ (* dy 2) (* y -1)))
                        1 1 0.001) 1000)
=> 2.716923932235896