SICP: Exercise 1.29

1.3 節「高階手続きによる抽象」から、ついに lambda が登場する。SICP を始める前には主に JavaScript を使っていたのだけれど、『JavaScript: The Good Parts』で匿名関数や高階関数の強力さが力説されているのを見てきたので、この辺りはかなり楽しく進めていた記憶がある。

というわけで、問題 1.29。

問題 1.29

Simpson の公式は上に示したのより更に正確な数値積分法である。Simpson の公式を使えば、\( a \) から \( b \) までの関数 \( f \) の積分は、偶整数 \( n \) に対し、\( h = \frac{b - a}{n} \) また \( y_k = f(a + kh) \) として $$ \frac{h}{3} [ y_0 + 4y_1 + 2y_2 + 4y_3 + 2y_4 + \cdots + 2y_{n - 2} + 4y_{n - 1} + y_n ] $$ で近似出来る。（\( n \) の増加は近似の精度を増加する。）引数として \( f \), \( a \), \( b \) と \( n \) をとり、Simpson 公式を使って計算した積分値を返す手続きを定義せよ。その手続きを使って（\( n = 100 \) と \( n = 1000 \) で）0 から 1 まで cube を積分し、またその結果を上の integral 手続きの結果と比較せよ。

「計算機プログラムの構造と解釈 第二版 1.3.1 引数としての手続き」より

高階手続きの導入に伴って、これまでのものより抽象的な手続きが出てきた。

sum は総和 $$ \sum_{n = a}^b\,f(n) = f(a) + \cdots + f(b) $$ を表す手続き。

;; term および next は手続き
(define (sum term a next b)
  (if (> a b)
      0
      (+ (term a)
         (sum term (next a) next b))))

inc は引数を 1 増やす手続き。

(define (inc n) (+ n 1))

恒等手続き。

(define (identity x) x)

定積分は細かく区間を分けて長方形で近似する。\( \Delta x \) を小さな実数として $$ \int_a^b f(x)\,dx \approx \left\{ f \left( a + \frac{\Delta x}{2} \right) + f \left( a + \Delta x + \frac{\Delta x}{2} \right) + f \left( a + 2\Delta x + \frac{\Delta x}{2} \right) + \cdots \right\} \Delta x $$

手続きとして書くと

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

解答

\( n \) は偶数のみを想定している。奇数のときにエラーとする手もありそうだけれど、係数が 4, 2, 4, … となるだけなので、チェックは省いた。

(define (integral-simpson f a b n)
  (let ((h (/ (- b a) n)))
    (define (y k) (f (+ a (* k h))))
    (define (term k)
      (* (cond ((or (= k 0) (= k n)) 1.0)
               ((even? k) 2.0)
               (else 4.0))
         (y k)))
    (* (/ h 3) (sum term 0 inc n))))

integral との比較

0 から 1 まで cube を積分した結果で、integral と integral-simpson を比較した。正確な値は 0.25。

;; 積分値、誤差を計算する

;; n = 100
(integral cube 0 1 0.01)
; 0.24998750000000042
(- (integral cube 0 1 0.01) 0.25)
; -1.249999999958229e-5

(integral-simpson cube 0 1 100)
; 0.24999999999999992
(- (integral-simpson cube 0 1 100) 0.25)
; -8.326672684688674e-17

;; n = 1000
(integral cube 0 1 0.001)
; 0.249999875000001
(- (integral cube 0 1 0.001) 0.25)
; -1.2499999899051595e-7

(integral-simpson cube 0 1 1000)
; 0.2500000000000002
(- (integral-simpson cube 0 1 1000) 0.25)
; 2.220446049250313e-16

長方形近似と比べると、Simpson の公式を使えばかなり精度が上がることが分かる。面白いのが、\( n = 100 \) のときの方が \( n = 1000 \) のときより真の値に近いこと。丸め誤差が蓄積したのだろうか。