[SICP] 問題 1.29 : Simpsonの公式

Simpsonの公式は上に示したのより更に正確な数値積分法である.
Simpsonの公式を使えば, aからbまでの関数fの積分は, 偶整数nに対し, h = (b - a)/nまたy_k = f(a + kh)として

(h/3)[y₀ + 4y₁ + 2y₂ + 4y₃ + 2y₄ + ・・・ + 2y_n-2 + 4y_n-1 + y_n]

で近似出来る. (nの増加は近似の精度を増加する.)
引数としてf, a, bとnをとり, Simpson公式を使って計算した積分値を返す手続きを定義せよ. その手続きを使って(n = 100とn = 1000で)0から1までcubeを積分し, またその結果を上のintegral手続きの結果と比較せよ.

Simpsonの公式の大カッコの中身をじっと見つめると, 次のようなパターンが見えてきます.

(y₀ + 4y₁ + y₂) + (y₂ + 4y₃ + y₄) + ・・・ + (y_n-4 + 4y_n-3 + y_n-2) + (y_n-2 + 4y_n-1 + y_n)

y₂に注目してください. 公式では係数が2となっていますが, それを2つに分けています.

つまり, 大カッコの中は３つの項をひとまとめにするとスッキリと書けます. 手続きを定義します.

(require (lib "racket/trace.ss"))

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

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

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

(define (simpson f a b n)
  (define h (/ (- b a) n))
  (define (y k) (f (+ a (* k h))))
  (define (term k) (+ (y k) (* 4.0 (y (+ k 1))) (y (+ k 2))))
  (define (next k) (+ k 2))
  (* (/ h 3.0)
     (sum term 0 next (- n 2))))

0から1までのcubeを積分してみます.

ようこそ DrRacket, バージョン 6.1 [3m].
言語: Pretty Big; memory limit: 2048 MB.
> (integral cube 0 1 0.01)
0.24998750000000042
> (integral cube 0 1 0.001)
0.249999875000001
> (simpson cube 0 1 100)
0.2500000000000001
> (simpson cube 0 1 1000)
0.2500000000000001
> 

計算結果を見ると精度が全く異なっています.