計算機プログラムの構造と解釈 第二版 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とでました。
