SICP 問題 1.29 （Simpsonの公式を使った積分計算）

Scheme SICP

【問題】

Simpsonの公式は上に示したのより更に正確な数値積分法である。
Simpsonの公式を使えば、a から b までの関数 f の積分は、



偶整数 n に対し、
h = (b-a)/n
また、


y_k = f(a+kh) として、
(h/3)[y_0 + 4・y_1 + 2・y_2 + 4・y_3 + 2・y_4 + ... + 2・y_(n-2) + 4・y_(n-1) + y_n]
で近似できる。（n の増加は近似の精度を増加する。）
引数として f, a, b, と n をとり、Simpson公式を使って計算した積分値を返す手続きを定義せよ。
その手続きを使って（n=100 と n=1000で）0 から 1 までcubeを積分し、またその結果を上のintegral手続きの結果と比較せよ。

【解答】

う〜ん、思考プロセスをログに取っていかなくても、すんなりできちゃったなぁ。。
とりあえずコード中にその時考えたことを記録しておこうか。

(define (simpson f a b n)
  ;yの関係式だけでも式にすると結構面倒なので内部手続き化してしまう。
  ;この時注意したのは、
  ;　「hの中で使用されているn」と「yの添字として渡されるk」は別物
  ;だということ。
  (define (y k)
    (f (+ a (* k (/ (- b a) n)))))

  ;名前付き再帰処理を定義。別に内部手続き化しても同じだからお好きな方で。
  (let loop ((k 0)       ;現在処理中の項のインデックスを表現。
	     (answer 0)) ;現在の加算結果（初期値は0とする）
    (cond ((= k 0)
	   ;初期値の場合は係数なしでyを呼ぶ。
	   (loop (+ k 1) (+ answer (y k)))) 
	  ((= k n)
	   ;最終値の計算の場合は係数なしでyを呼び、加算後にh/3を乗算する。
	   (* (/ (/ (- b a) n) 3) (+ answer (y k))))
	  ((odd? k)
	   ;奇数の場合は係数を4にしてyを呼んで加算。
	   (loop (+ k 1) (+ answer (* 4 (y k)))))
	  (else
	   ;偶数の場合は係数を2にしてyを呼んで加算。
	   (loop (+ k 1) (+ answer (* 2 (y k))))))))

定義はこんな感じ。んで、cubeをfとして使いたいのでこれも定義してしまおう。

(define (cube x)
  (* x x x))

さて、この２つの手続きを読み込ませて実行してみる！



gosh> (simpson cube 0 1 100)

1/4

gosh> (simpson cube 0 1 1000)

1/4

gosh>
お〜、ドンピシャ。Gaucheの正確数だっけ？有理数で結果が取得できるぞ。

問題文中には



integral手続きの結果と比較せよ。
とあるので、この問題の前に定義されていたintegral手続きとそれから呼び出されるされるsum手続きもここに記載する。

(define (integral f a b dx)
  (define (add-dx x)
    (+ x dx))
  (* (sum f (+ a (/ dx 2.0)) add-dx b)
     dx))

(define (sum term a next b)
  (if (> a b)
      0
      (+ (term a)
	 (sum term (next a) next b))))

こいつらを読み込ませて、実行してみる。
ちなみにこいつらは小さい値(0.01とか0.001とか）を計算の元ネタとして積分していく。



gosh> (integral cube 0 1 0.01)

0.24998750000000042

gosh> (integral cube 0 1 0.001)

0.249999875000001

gosh> (integral cube 0 1 (/ 1 1000))

0.249999875000001

gosh>
というわけで、その小さい値を有理数にしても、得られる結果は非正確数（だっけ？）になる。
なので、Simpsonの公式を使った結果の方が正確な値を取得できる。