プログラムで定積分するだけ

16403 ワード
積分 dlang 数学 dlang テキストリンク

どうでもいい動機

宿題が面倒だった。

数学Ⅱ(高校)の最後を飾る分野は、単純計算が面倒くさい微分・積分法である。特に定積分では分数の嵐に襲われ、ひたすら分数の足し算をしていた小学校の記憶が蘇るほど面倒くさい。頭をちょっとでもひねる問題ならともかく、問題集の最初の方なんてやってられるか！と思って書いてみた。
でも結局宿題は普通にやった。

積分について

一応おさらい

数学Ⅱでは、定積分をこのように定義している。

関数f(x)の原始関数の一つをF(x)としたとき、
\\ \int_a^b f(x) dx=\left[F(x)\right]_a^b=F(b)-F(a)

幾何学的には、これは

グラフy=f(x),y=0,x=a,x=bによって囲まれた領域の面積は
\\\int_a^b f(x) dx \\(ただしa\leq x \leq bにおいてf(x)\geq 0)

ということにもなる。

アプローチ

上記の幾何学的な云々は、「aからbの間に微小区間を取り($\Delta x$とか)、それと$f(\Delta x)$を掛け合わせてできた細い短冊を積み重ねていく」という解釈で説明されることが多い。

ちょうどいい図解があったので引用
http://oto-suu.seesaa.net/article/459128337.html

とりあえずこれを使って実装してみようと思う。

今回は$f(x)=3x^2-6x+5$ を $0\leq x\leq3$で積分する。
値は、$\int_0^3 f(x) dx=\left[ x^3-3x^2+5x\right]_0^3=15$

実装

とりあえず書いてみる

やっぱりD言語がナンバーワン！

ver1.d
import std.stdio;
import std.string;
import std.conv;

void main(){
    auto f=(real x)=>3*x^^2-6*x+5;

    real trueVal=15;
    real bottom=0;
    real top=3;

    "enter times to split:".write;
    auto spl=readln.chomp.to!int;
    foreach(s;1..spl+1){
        real delta=(top-bottom)/s;
        real sum=0;
        foreach(n;1..s+1){
            sum+=f(bottom+delta*n)*delta;
        }
        writef("true:%s | delta:%12s | approximate:%12s | error:%s\n",
                trueVal,delta,sum,sum/trueVal);
    }   
}

最初に区間を何分割するか(Δxをどのくらい細かくするか)を入力し、1からその回数まで1ずつ増やして計算していく。
出力のapproximateが近似値、errorは真値との誤差(今回は1に近いほうがいい)。

100割までしてみた。

実行結果
enter times to split:100
true:15 | delta:           3 | approximate:          42 | error:2.8
true:15 | delta:         1.5 | approximate:      25.125 | error:1.675
true:15 | delta:           1 | approximate:          21 | error:1.4
true:15 | delta:        0.75 | approximate:     19.2188 | error:1.28125
true:15 | delta:         0.6 | approximate:       18.24 | error:1.216
#中略
true:15 | delta:   0.0309278 | approximate:     15.1406 | error:1.00937
true:15 | delta:   0.0306122 | approximate:     15.1392 | error:1.00928
true:15 | delta:    0.030303 | approximate:     15.1377 | error:1.00918
true:15 | delta:        0.03 | approximate:     15.1363 | error:1.00909

まあそこそこ。
ただこの例だと

f(a+\Delta x)\Delta x+f(a+2\Delta x)\Delta x+...+f(a+s\Delta x)\Delta x
\\=\{f(a+\Delta x)+f(a+2\Delta x)+...+f(a+s\Delta x)\}\Delta x

左辺のような計算をしていて非効率なので右辺に書き換えよう。

ver1.d
//省略
foreach(n;1..s+1){
    sum+=f(bottom+delta*n);
}
sum*=delta;
//省略
実行結果
#中略
true:15 | delta:    0.030303 | approximate:     15.1377 | error:1.00918
true:15 | delta:        0.03 | approximate:     15.1363 | error:1.00909

結果は変わらないことを確認。よしよし。
ちなみに10万回やったら誤差1.00001くらいになった。さすがreal型、80bit精度は一味違いますな(知らんけど)。でもそれだとこの後改善しても解りづらいので、とりあえずこれを基準とする。

この例だと、さっきの図でいうΔxの右端部分でf(x)をとっている。

中央で取ったほうが誤差は少ないかもしれない。

真ん中でとる

一箇所変えるだけ

ver2.d
//省略
sum+=f(bottom+delta*(n-0.5));
//省略

これで$\Delta x$の真ん中がとれるはず。

実行結果
enter times to split:100
true:15 | delta:           3 | approximate:        8.25 | error:0.55
true:15 | delta:         1.5 | approximate:     13.3125 | error:0.8875
true:15 | delta:           1 | approximate:       14.25 | error:0.95
true:15 | delta:        0.75 | approximate:     14.5781 | error:0.971875
true:15 | delta:         0.6 | approximate:       14.73 | error:0.982
#中略
true:15 | delta:   0.0309278 | approximate:     14.9993 | error:0.999952
true:15 | delta:   0.0306122 | approximate:     14.9993 | error:0.999953
true:15 | delta:    0.030303 | approximate:     14.9993 | error:0.999954
true:15 | delta:        0.03 | approximate:     14.9993 | error:0.999955

絶対値で比較すると、有効数字を揃えても前回から誤差が0.0905縮まった。10%と考えると結構大きい。

台形近似

世の中には台形近似というものもあるらしい。簡単に言うと先程の短冊が台形になったものだ。当然台形の面積は底辺×高さ÷2で求められるので、これを使った方が形的には近そう。
計算の処理自体は

\{\frac{f(a)+f(a+\Delta x)}{2}+\frac{f(a+\Delta x)+f(a+2\Delta x)}{2}+...+\frac{f(a+(s-1)\Delta x)+f(a+s\Delta x)}{2}\}\Delta x\\
=\{\frac{f(a)}{2}+f(a+\Delta x)+f(a+2\Delta x)...f(a+(s-1)\Delta x)+\frac{f(a+s\Delta x)}{2}\}\Delta x

とまとめられる。

ver3.d
import std.stdio;
import std.string;
import std.conv;

void main(){
    auto f=(real x)=>3*x^^2-6*x+5;

    //integral f(x) 0->3
    real trueVal=15;
    real bottom=0;
    real top=3;

    "enter times to split:".write;
    auto spl=readln.chomp.to!int;
    foreach(s;1..spl+1){
        real delta=(top-bottom)/s;
        real sum=0;

        sum+=f(bottom)/2;
        foreach(n;1..s){
            sum+=f(bottom+delta*n);
        }
        sum+=f(bottom+delta*s)/2;

        sum*=delta;
        writef("true:%s | delta:%12s | approximate:%12s | error:%s\n",
                trueVal,delta,sum,sum/trueVal);
    }   
}
実行結果
enter times to split:100
true:15 | delta:           3 | approximate:        28.5 | error:1.9
true:15 | delta:         1.5 | approximate:      18.375 | error:1.225
true:15 | delta:           1 | approximate:        16.5 | error:1.1
true:15 | delta:        0.75 | approximate:     15.8438 | error:1.05625
true:15 | delta:         0.6 | approximate:       15.54 | error:1.036
#中略
true:15 | delta:   0.0309278 | approximate:     15.0014 | error:1.0001
true:15 | delta:   0.0306122 | approximate:     15.0014 | error:1.00009
true:15 | delta:    0.030303 | approximate:     15.0014 | error:1.00009
true:15 | delta:        0.03 | approximate:     15.0013 | error:1.00009

...あれ、精度下がった？
前回が0.999955、今回が1.00009だから、桁を揃えても0.0003の誤差増加
なぜだろうか。

まとめ

なんか知らんが台形近似だと誤差が下がった。しょうもないコーディングミスかもしれないし、もっと深いわけがあるのかもしれない。
しかしとりあえず数Bの宿題がたくさんあるので、それが終わったら調べてみる。