SICP 3.5.4 - 3.5.5

3.77

#lang racket
(define (integral delayed-integrand initial-value dt)
  (let ((integrand (froce delayed-integrand)))
    (cons-stream initial-value
      (if (stream-null? integrand)
          the-empty-stream
          (integral (stream-cdr integrand)
                    (+ (* dt (stream-car integrand))
                       initial-value)
                    dt)))))

3.78

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

3.79

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

3.80

(define (RLC R L C dt)
  (define (rlc vC0 iL0)
    (define vC (integral (delay dvC) vC0 dt))
    (define iL (integral (delay diL) iL0 dt))
    (define dvC (scale-stream iL (/ -1 C)))
    (define diL (add-streams
                  (scale-stream vC (/ 1 L))
                  (scale-stream iL (- (/ R L)))))
    (stream-map (lambda (x y) (cons x y)) vC iL))
  rlc)

(define RLC1 (RLC 1 1 0.2 0.1))

(display-stream-head (RLC1 10 0) 10)

結果

(10 . 0)
(10 . 1.0)
(9.5 . 1.9)
(8.55 . 2.66)
(7.220000000000001 . 3.249)
(5.5955 . 3.6461)
(3.77245 . 3.84104)
(1.8519299999999999 . 3.834181)
(-0.0651605000000004 . 3.6359559)
(-1.8831384500000004 . 3.2658442599999997)
'done

実行には以下のhelperを使う

#lang racket
(require (prefix-in strm: racket/stream))

;helper
(define-syntax cons-stream
  (syntax-rules ()
    ((_ a b) (strm:stream-cons a b))))
(define stream-car strm:stream-first)
(define stream-cdr strm:stream-rest)
(define stream-null? strm:stream-empty?)
(define (scale-stream stream factor)
  (stream-map (lambda (x) (* x factor)) stream))
(define the-empty-stream strm:empty-stream)
(define (display-stream s) (stream-for-each display-line s))
(define (add-streams s1 s2) (stream-map + s1 s2))
(define (stream-map proc . argstreams)
  (if (stream-null? (car argstreams))
      the-empty-stream
      (cons-stream
       (apply proc (map stream-car argstreams))
       (apply stream-map
              (cons proc (map stream-cdr argstreams))))))
(define (display-stream-head s n)
  (define (iter s n)
    (if (<= n 0)
      'done
      (begin
        (display (stream-car s))
        (newline)
        (iter (stream-cdr s) (- n 1)))))
  (iter s n))
(define (display-line x)
  (newline)
  (display x))

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

3.81

入力として ‘generate が来たら乱数を生成。

数値が来たらその数値を使ってresetという動作。

3.6では代入で値を保持していたが、ストリームの中に閉じるので代入は不要になる。

ストリームの中で値を都度 last-value として次の処理に渡してやれば良い。

(define (rand-stream s seed)
  (define (generate last-val)
    ;lcdとか使うと乱数っぽくなるが、テストを楽にするために簡略化
    (+ last-val 1))
  (define (reset new-seed) new-seed)
  (define (dispatch req last-val)
    (cond ((eq? req 'generate) (generate last-val))
          ((number? req) (reset req))
          (else (error "unknown message" req))))
  (define (iter req-stream last-val)
    (let ((new-val (dispatch (car req-stream) last-val)))
      (cons-stream
        new-val
        (iter (stream-cdr req-stream) new-val))))
  (iter s seed))

;test
(define test-stream (list 'generate 5 'generate 'generate 10 20 'generate))

(display-stream-head (rand-stream test-stream 1) 7)

結果。（実行時は3.80と同じhelperを追加する）

2
5
6
7
10
20
21
'done

3.82

試行毎に乱数生成していたのを、あらかじめ乱数ストリームを渡すように改造する。

(define (estimate-integral p x1 x2 y1 y2 )
  (define x-stream (random-in-range-stream x1 x2))
  (define y-stream (random-in-range-stream y1 y2))
  (stream-map
   (lambda (p) (* (* (- x2 x1) (- y2 y1))p))
   (monte-carlo (stream-map p x-stream y-stream) 0 0)))


;単位円をテストする関数
(define (circle x y)
   (<= (+ (square x) (square y)) 1.0))

;中心 (50,50) 半径 50 の円
(define (circle2 x y)
  (<= (+ (square (- x 50.0)) (square (- y 50.0))) (square 50.0)))

; piの見積もり (pi=S/r^2)
; S 円の面積
; r 半径
(define (estimate-pi-stream S-stream r)
  (stream-map (lambda (s) (/ s (square r))) S-stream))

;test
(display-stream-head
 (estimate-pi-stream (estimate-integral circle -1 1 -1 1) 1.0) 100)

(display-stream-head
 (estimate-pi-stream (estimate-integral circle2 0 100 0 100) 50.0) 100)

実行結果（長いので最後の要素だけ）

;単位円
...
2.72
'done
;半径50の円
3.16
'done

精度がいまいちなのは乱数の質の問題？

なお、実行には以下のhelperを使う

#lang racket
(require (prefix-in strm: racket/stream))

;helper
(define-syntax cons-stream
  (syntax-rules ()
    ((_ a b) (strm:stream-cons a b))))
(define stream-car strm:stream-first)
(define stream-cdr strm:stream-rest)
(define stream-null? strm:stream-empty?)
(define (scale-stream stream factor)
  (stream-map (lambda (x) (* x factor)) stream))
(define the-empty-stream strm:empty-stream)
(define (display-stream s) (stream-for-each display-line s))
(define (add-streams s1 s2) (stream-map + s1 s2))
(define (stream-map proc . argstreams)
  (if (stream-null? (car argstreams))
      the-empty-stream
      (cons-stream
       (apply proc (map stream-car argstreams))
       (apply stream-map
              (cons proc (map stream-cdr argstreams))))))
(define (display-stream-head s n)
  (define (iter s n)
    (if (<= n 0)
      'done
      (begin
        (display (stream-car s))
        (newline)
        (iter (stream-cdr s) (- n 1)))))
  (iter s n))
(define (display-line x)
  (newline)
  (display x))
(define ones (cons-stream 1.0 ones))
(define (square x) (* x x))
(define (random-in-range low high)
  (let ((range (- high low)))
    (+ low (random range))))
(define (random-in-range-stream low high)
  (stream-map
   (lambda (x) (random-in-range low high))
   ones))

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