SICPを読む(20) 1.3.1(2) 区分求積法による積分

さて、昨日のアレだ。アレ。

$\displaystyle \int_{a}^{b}f = \left[f(a + \frac{dx}{2})+f(a + dx + \frac{dx}{2})+f(a + 2dx + \frac{dx}{2})+\cdots \right]dx$

なんか難しく書いてあるけど、簡単だった。

積分は面積だ。

まず、 $x^2$ のグラフを描こう。

面積を計れ。

・・・。

       *
       *
      **
      **
     ***
   *****
********
+------+
0      1

紙面の関係上 $x^2$ にはなってない。

"分割してたしざん"これこそが積分だった。

              *
              *
            * *
            * * = x^2の面積。
          * * *
      * * * * *
* * * * * * * *

もう一度式を見たら、意味不明な式が面積に見えてくるはず。

"分割してたしざん"

が見えればOK。

で、一般的な式に直すと、

$\Delta x = \frac{b - a}{n}$
$x_k = (a + k\Delta x)$
$\int_{a}^{b}f(x)dx = \lim_{n\to\infty}\sum_{k=0}^{n-1}f(x_k)\Delta x$

うは。なんか凄い式。

$f(x_k)\Delta x$ で、1個の面積。全部足せば、面積が出る。

これが区分求積法。太古の昔から使われていたらしい。

面倒なのは、SICPの求め方。

$x_k = (a + k\Delta x + \frac{\Delta x}{2})$
$\int_{a}^{b}f(x)dx = \lim_{n\to\infty}\sum_{k=0}^{n-1}f(x_k)\Delta x$

元に戻ると、

$\displaystyle \int_{a}^{b}f = \left[f(a + \frac{dx}{2})+f(a + dx + \frac{dx}{2})+f(a + 2dx + \frac{dx}{2})+\cdots \right]dx$

こんな感じ。

"分割してたしざん"

これが積分。

僕は、記号の山を見てびびってた。

「記号は抽象」ってSICPに書いてあった。

抽象の壁を破ったら、積分なんてメチャクチャ簡単だった!!

では、Schemeにする。

(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) (+ dx x))
  (* (sum f (+ a (/ dx 2.0)) add-dx b)
     dx))

(define (cube n) (* n n n))
(integral cube 0 1 0.01) ; 0.24998

プログラム見ると安心するね。要は足し算と掛算。

さて、積分を抽象化しました。

凄いのはここから。

(integral sqrt 0 1 0.01) ; 0.6667
(integral exp  0 1 0.01) ; 1.7182
(integral sin  0 3.14 0.01) ; 2.0000
(integral asin 0 1 0.01) ; 0.5707

積分の抽象化によってあらゆる関数の積分が出来るようになる。

抽象が抽象を生み、その抽象が抽象を生む。

関数型言語を学ぶ理由はここにあると思っている。やはりSICPを読まなければならない。

問題は数学・・・。

追記

nの計算を勘違いしてた。

|(sum #<procedure:cube> 0.05 #<procedure:add-dx> 1)
| (sum #<procedure:cube> 0.15000000000000002 #<procedure:add-dx> 1)
(略
| | | | | (sum #<procedure:cube> 0.9499999999999998 #<procedure:add-dx> 1)
| | | |[10](sum #<procedure:cube> 1.0499999999999998 #<procedure:add-dx> 1)
| | | |[10]0
| | | | | 0.8573749999999997
(略
|2.4874999999999994
0.24874999999999994

(+ a (/ dx 2.0))してるので、dxが0.1の時、終端は約1.05になり、確実に1を越える。

aのみにすると、誤差があるので、1を越えたり、越えなかったり、意味不明な値になる。

| | | |[10](sum #<procedure:cube> 0.9999999999999999 #<procedure:add-dx> 1) ← 10分割
| | | |[11](sum #<procedure:cube> 1.0000000000000002 #<procedure:add-dx> 1) ← 11分割

超危険。終端bに補正値を入れる必要がある。

ふつうの区分求積法bに(/ dx 2)の補正値を追加する。

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

確実に1を計算しなくなる。n-1を実現できた。

精度はかなり悪い。