計算機プログラムの構造と解釈 第二版 P33 問題1.29
はいはい、そうそう、
ちょっと前までの説明の
dx <-(増分のこと)
の代わりに、
h = (b - a)/n
ととかやって、
y[k] = f(a + k * h)
とかやる訳ね。
よゆうよゆう
と余裕をかましていたのでした。
数分後 - - -
ああ、もうわかる訳がない!!
こんな問題!!!
なんだこりゃ、、、、
h/3(y[0] + 4y[1] + 2y[2] + 4y[3] + 2y[4] ... + 2y[n-2] + 4y[n-1] + y[n])
いや、はじめはね、
nが奇数の時、4掛けて、
nが偶数の時、2掛けりゃいいと思った。
だけど、、おいおい、なんだよ、初項と終項、どうなってんだよ!!!
と、怒りに燃えたのでした。
もうしょうがないのでここを解決するために、
せっかくここまでがんばったのですが、、、、、、グーグル先生に聞いたのでした。
断腸の思いだ!!!。
素人くさいSICP読書会 Wiki
http://www.csus4.net/hiki/SICPReading/?1.29
公式中の項2つで term を定義し条件分岐を排除
;;( (2*y[0] + 4*y[1]) + ... + (2*y[n-2] + 4*y[n-1])) + 2*y[n] - y[0] - y[n]
すごい、、なるほど!!
頭の良い人というのはいるものです。
説明いりますか?(未来の自分に向けて)
あ、説明?必要?ですねー。(どうせ明日の僕は忘れているので、明日見てもわかるように。)
どういうコトかというと、とりあえずはこうしちゃう。
(2y[0] + 4y[1] + 2y[2] + 4y[3] + 2y[4] ... + 2y[n-2] + 4y[n-1] + 2y[n])
求めたいのは
(y[0] + 4y[1] + 2y[2] + 4y[3] + 2y[4] ... + 2y[n-2] + 4y[n-1] + y[n])
だから、引き算をしてみましょう。
(2y[0] + 4y[1] + 2y[2] + 4y[3] + 2y[4] ... + 2y[n-2] + 4y[n-1] + 2y[n]) - ( y[0] + 4y[1] + 2y[2] + 4y[3] + 2y[4] ... + 2y[n-2] + 4y[n-1] + y[n]) --------------------------------------------------------------------------- y[0] + y[n]
だので、余分な部分を引いてしまえば同じ訳です。
( y[0] + 4y[1] + 2y[2] + 4y[3] + 2y[4] ... + 2y[n-2] + 4y[n-1] + y[n]) <-求めたい部分 = (2y[0] + 4y[1] + 2y[2] + 4y[3] + 2y[4] ... + 2y[n-2] + 4y[n-1] + 2y[n]) - y[0] - y[n]
んで、次に思うのは、奇数の時と偶数の時を考えるんだけど、そうすると、
if odd? n
みたいな感じで書かないといけない。
それはいやなので、nを2ずつインクリメントさせていけばいい
そうすると、
0の時は、
( n (つまり第0項)の結果 ) * 2
( n + 1 (つまり第1項)の結果 ) * 4
インクリメントして2の時は、
( n (つまり第2項)の結果 ) * 2
( n + 1 (つまり第3項)の結果 ) * 4
って感じでうまくいきそう。
でもそうすると、最後の第n項はうまくいかないんだよね。
何でかっていうと、sumの計算の終了条件が
(> a b)
だからなんだけど、
たとえば、nを1000までインクリメントさせていくと、
nが1000になったとき、
(> 1000 1000)
では終了しできない。(sumは書き換えないのが前提)
そうすると、2ずつインクリメントされているから、
( n (つまり第1000項)の結果 ) * 4 ( n + 1 (つまり第1001項)の結果 ) * 2
まで計算されて、
nが1002の場合までsumをとられて終了する。
それじゃまずい。
結果として、この場合、
終了条件に合うように、あらかじめnで表される終項を
(- n 2)
としておいてやるのだ。
そうすると、終了条件が998ということになるので、
結果、998の時は
(> 998 998)
は通って
( n (つまり998項目)の結果 ) * 4 ( n + 1 (つまり999項目)の結果 ) * 2
までは計算されて、
nが1000のときsumの終了条件
(> 1000 998)
に引っかかって止まる。
そうすると998、999項は計算されるが
第1000項は計算されないので、後でたしてやればよい。
ということで、
前述の
;;( (2*y[0] + 4*y[1]) + ... + (2*y[n-2] + 4*y[n-1])) + 2*y[n] - y[0] - y[n]
が、きれいに理解できる。
ので結局
(((2*y[0] + 4*y[1]) + ... + (2*y[n-2] + 4*y[n-1])) + y[n]) - y[0]
みたいな。
ここまでわかればあとは、余裕のよっちゃんだね。
ということで、黙々と作業を続けできたのがこれだ!!
#!/usr/local/bin/gosh ;; -*- coding: utf-8 -*- (use ggc.debug.trace) (use math.mt-random) (define (sum term a next b) (if (> a b) 0 (+ (term a) (sum term (next a) next b)))) (define (cube a) (* a a a)) ;;integral (define (integral f a b dx) (define (add-dx x) (+ x dx)) (* (sum f (+ a (/ dx 2.0)) add-dx b) dx)) ;;Simpson (define (simpson f a b n) (define h (/ (- b a) n)) (define (yk k) (f (+ (* k h) a))) (define (next n) (+ n 2)) (define (term k) (+ (* 2.0 (yk k)) (* 4.0 (yk (+ k 1))))) (/ (* h (- (+ (sum term 0 next (- n 2)) (yk n)) (yk 0))) 3.0)) ;; main (define (main args) (display "simpsonを利用して0から1までのcubeを、100個に区切って計算 : ") (display (simpson cube 0 1 100))(newline) (display "simpsonを利用して0から1までのcubeを、1000個に区切って計算 : ") (display (simpson cube 0 1 1000))(newline)(newline) (display "integralを利用して0から1までのcubeを、100個(0.01刻み)に区切って計算 : ") (display (integral cube 0 1 0.01))(newline) (display "integralを利用して0から1までのcubeを、1000個(0.001刻み)に区切って計算 : ") (display (integral cube 0 1 0.001))(newline) 0)
integralとの結果比較(実行結果)
simpsonを利用して0から1までのcubeを、100個に区切って計算 : 0.25000000000000006 simpsonを利用して0から1までのcubeを、1000個に区切って計算 : 0.25 0と1の間のcubeの積分(0.01刻み): 0.24998750000000042 0と1の間のcubeの積分(0.001刻み): 0.249999875000001
Simpsonの方が本当に精度が高いね。
1000だと、もう誤差が小さすぎて、0.25とでました。