問題3.77-3.82

ストリームの続きです。
フィードバックループのあるシステムのストリームによるモデル化では、cons-streamの中に隠れたdelayだけでは表現できないので、別に明示的にdelayを使う必要がある、ということらしいです。たぶん。

3.77

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

integralの第一引数を遅延オブジェクトにし、streamの二番目以降でforceして使うようにする。integralを再帰的に呼び出すときに、integrandのstream-cdrをdelayするのを忘れると無限ループになるので注意。
(一回忘れました。)

3.78

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

たぶんこんな感じ。y(t)=e^tとして試してみます。

> (define s (solve-2nd 0 1 0.0001 1 1))
> (stream-ref s 10000)
2.7181459268252266
> (define s (solve-2nd 1 0 0.0001 1 1))
> (stream-ref s 10000)
2.7181459268252266

で正しいように思います。

3.79

(define (solve-2 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)

例のごとくe^tで確認します。

> (define s (solve-2 (lambda (x y) (+ (* 0 x) (* 1 y))) 1 1 0.001))
> (stream-ref s 1000)
2.716923932235896

おっけーのようです。簡単にかけていいですね。

3.80

対のストリームという文をどっちに解釈したらいいのか分からないんだけど、こっちでいいのかな。

(define (RLC r l c dt)
  (lambda (v_c0 i_l0)
    (define v_c (integral (delay dv_c) v_c0 dt))
    (define i_l (integral (delay di_l) i_l0 dt))
    (define dv_c (scale-stream i_l (- (/ 1 c))))
    (define di_l (add-streams (scale-stream v_c (/ 1 l))
                              (scale-stream i_l (- (/ r l)))))
    (cons v_c i_l)))

時間経過による回路の状態を追ってみます。

> (define rlc1 ((RLC 1 1 0.2 0.1) 10 0))
> (stream-ref (car rlc1) 10)
-3.5160605800000004
> (stream-ref (cdr rlc1) 10)
2.750945989
> (stream-ref (car rlc1) 1000)
4.5996318176662634e-11
> (stream-ref (cdr rlc1) 1000)
-2.1020456106564936e-11
> (stream-ref (car rlc1) 100000)
4.9406564584125e-324
> (stream-ref (cdr rlc1) 100000)
-4.9406564584125e-324

無限時間経過でv_cもi_lも0に近づくようです。
100くらい連続で値を出してみると、どうやら減衰振動しているみたいです。
回路の知識に乏しいのであんまり自信ないですが、キャパシタの両極の電荷が行ったりきたりしながら同じになっていくイメージで考えてみるとまあ合ってる、ような気がします。

3.81

問題文の意味が分からないので省略

3.82

(define (rand-update x) (modulo (+ (* 1000 x) 100000) 104729))
(define random-init 8)
(define random-numbers
  (cons-stream random-init
               (stream-map rand-update random-numbers)))
(define random-init2 123)
(define random-numbers2
  (cons-stream random-init2
               (stream-map rand-update random-numbers2)))
(define (monte-carlo experiment-stream passed failed)
  (define (next passed failed)
    (cons-stream
     (/ passed (+ passed failed))
     (monte-carlo
      (stream-cdr experiment-stream) passed failed)))
  (if (stream-car experiment-stream)
      (next (+ passed 1) failed)
      (next passed (+ failed 1))))
(define (estimate-integral P x1 x2 y1 y2)
  (define (random-in-range-s low high s)
    (let ((range (- high low)))
      (stream-map (lambda (x) (+ low (modulo x range))) s)))
  (define width-random (random-in-range-s x1 x2 random-numbers))
  (define hight-random (random-in-range-s y1 y2 random-numbers2))
  (define test-stream
    (stream-map P width-random hight-random))
  (define (rect x1 x2 y1 y2) (* (- x2 x1) (- y2 y1)))
  (stream-map (lambda (s) (* (rect x1 x2 y1 y2) s)) (monte-carlo test-stream 0 0)))
(define estimate-pi
  (stream-map (lambda (s) (/ s (sqr 2)))
              (estimate-integral
               (lambda (x y) (<= (+ (sqr (- x 2)) (sqr (- y 2))) (sqr 2))) 0.0 10.0 0.0 10.0)))

なんというか、random回りがひどくやっつけ仕事ですが、一応ちゃんと動きます。

> (define s (estimate-integral 
             (lambda (x y) (<= (+ (sqr (- x 5)) (sqr (- y 7))) (sqr 3))) 2.0 8.0 4.0 10.0))
> (stream-ref s 1000)
26.613386613386613
> (stream-ref s 10000)
27.076492350764923

> (stream-ref estimate-pi 1000)
3.046953046953047
> (stream-ref estimate-pi 10000)
3.2496750324967505

やっぱり精度はよくないですね。実装が悪いのかモンテカルロ積分がもともとそういうものなのかはよく分かりませんが。


やっと3章が終わりました。予定の2倍くらい時間かかってますが。
きりがいいので今日はここまで。
次からは4章の超言語的抽象に入ります。