簡単な数値積分をプログラムする

カリカリ…、カチカチ…。なにやら家の中で音が聞こえます。ミケのようです。

ミケ：ん～、ここをこうすればいいのか。

どうやら今日はマロが散歩に行っているようで、ミケは一人みたいです。

ミケ：よ～し、できたぞぉ。なかなか良い近似だ。

ミケ：マロがいなくてさみしいけど、まぁたまにはこういうのもいいね。

ミケ：それじゃあ今回のブログには数値積分のプログラムについて書こう～

今回は数値積分のプログラムについてのようです。

予備知識

数値積分のプログラムの前に、数学についての知識を軽く説明していくね。

微分について

積分について説明する前に微分について説明しておこう。

微分っていうのは、ある関数の変化率を求める操作のことだよ。「微かに分かる」って意味じゃないから注意してね。

具体的には次の計算をするよ。

$${\lim_{h\to 0}\frac{f(x+h)-f(x)}{h}}$$

この極限を$${f'(x)}$$で表して、$${f(x)}$$の導関数というよ。このとき$${f'(a)}$$の値は$${f(a)}$$の値の傾きを表しているよ。

つまり導関数は元の関数の各点の傾きを集めた関数って考えるのがいいかもしれない。

例えば、$${f(x)=x^2}$$の導関数$${f'(x)}$$は

$${\lim_{h\to0}\frac{(x+h)^2-x^2}{h}=\lim_{h\to0}\frac{2xh+h^2}{h}=\lim_{h\to0}2x=2x}$$

で$${2x}$$となることがわかるよね。

もっとわかりやすい例をだすなら、$${f(x)=ax+b}$$の導関数を求めてみるといいかもしれないね。これは中学校で習う一次関数の一般式だから、傾きは常に$${a}$$になるはずだけど、果たして…?

計算してみると

$${\lim_{h\to0}\frac{a(x+h)+b-(ax+b)}{h}=\lim_{h\to0}\frac{ah}{h}=a}$$

この通り、$${a}$$になるね。

不定積分

$${a}$$に$${b}$$を足すと、$${a+b}$$になるよね。反対に、$${a+b}$$から$${b}$$を引くと$${a}$$に戻るよね？

このように互いが逆の演算になっているものを逆演算というよ。

足し算の逆演算は引き算、掛け算の逆演算は割り算だね。

じゃあ微分の逆演算って…?

結論から言うと、微分の逆演算は不定積分という演算だよ。

式で関係を表すと、

$${f(x)=\frac{d}{dx}\int f(x)dx}$$

となるよ。ちなみに$${\frac{d}{dx}f(x)=f'(x)}$$だよ。単に記号を変えたってだけさ。

$${F(x)=\int f(x)dx}$$とするとき、$${F(x)}$$は$${f(x)}$$の原始関数というよ。

例えば、$${f(x)=x}$$の原始関数を求めてみると、

$${\frac{d}{dx}x^2=2x}$$より、不定積分は微分の逆演算であることを利用して、

$${\int xdx=\frac{1}{2}\int 2xdx=\frac{1}{2}x^2+C}$$

となるよ。いきなり出てきたこの$${+C}$$というのは積分定数というもので、不定積分を計算する際は必ずつけないといけないよ。なぜなら、原始関数は定数の分だけ不定だからだよ(詳しくは割愛)。

定積分

ある関数$${y=f(x)\quad (a\leq x \leq b \to f(x)>0)}$$のx=aからbまでのグラフ下の領域の面積を求めるとしよう。

この時の面積を、

$${\displaystyle \int_a^b f(x)dx}$$

と表すことにする。これを定積分というよ。

また、

$${\displaystyle \int_a^b f(x)dx=\lim_{b\to\infty}\sum_{k=1}^n f(\zeta_k)\Delta x}$$

と表し、$${\displaystyle \sum_{k=1}^n f(\zeta_k)\Delta x}$$

を積和というよ。

とはいえ、いちいちこの公式を使って計算していったら面倒だよね？

そこで登場するのが微積分学の基本定理だ！

$${\displaystyle f(x)=\frac{d}{dx}\int_a^x f(t)dt}$$

この公式は定積分と不定積分の関係を表す式で、この式を用いることで、

$${\displaystyle \int_a^b f(x)dx=[F(x)]_a^b=F(b)-F(a)}$$

と定積分の値が、原始関数の積分上限から積分下限までを引いた値であることがわかるよね。

これを用いて試しに$${\int_0^{\pi}\sin x dx}$$を計算してみると、

$${\displaystyle \int\sin xdx=-\cos x+C}$$より、

$${\displaystyle \int_0^{\pi}\sin x dx=[-\cos x]_0^{\pi}=-\cos \pi +\cos 0=1+1=2}$$

であり、あの曲線の面積が2であることがわかるよ！

数値積分

数値積分っていうのは平たく言えば、「定積分の値を近似する」方法だよ！

原始関数が求められなくても、定積分の値を求めたい！ってときによく使うよ！

具体的には積和をどのように近似していくかで方法が変わるよ。詳しい説明は省いて今回は、公式だけ紹介するね。

台形公式

$${\displaystyle \int_a^b f(x)dx\approx \frac{h}{2}\{f(a)+2f(a+h)+2f(a+2h)+\cdots+2f(a+(n-1)h)+f(b)\}}$$

シンプソンの公式

$${\displaystyle \int_a^b f(x)dx\approx \frac{h}{3}\{f(a)+4f(a+h)+2f(a+2h)+\cdots+2f(a+(n-2)h)+4f(a+(n-1)h)+f(b)\}}$$

数値積分のプログラム

よし、じゃあ予備知識を頭に入れたところで、肝心のプログラムを見ていこう！今回もJavaを使ったよ！

3つのクラスで構成されているよ。

package javas.lab;

import javas.maths.*;

public class Main {
    public static void main(String[] args) {
        System.out.println(integral.simpson_rule(new Function() {
            @Override
            public double func(double x) {
                return 1/x;
            }   
            @Override
            public boolean isEvenFunction() {
                return false;
            }
            @Override
            public boolean isOddFunction() {
                return true;
            }
        }, 2, 6, 4));
    }

}

package javas.maths;

public interface Function {
    public double func(double x);
    // 偶関数かどうか
    public boolean isEvenFunction();
    // 奇関数かどうか
    public boolean isOddFunction();
}

package javas.maths;

public class Integral {
    /**
     * 
     * @param fx    :被積分関数
     * @param begin :積分下限
     * @param end   :積分上限
     * @param n :帯の数　多ければ多いほどよい近似値が得られる
     * @return 台形公式による定積分の近似値
     */
    public static double trapezoidal_rule(Function fx,double begin,double end,int n){
        if(begin==end)return 0;
        if(fx.isOddFunction()&&begin==-end)return 0;

        double sum=0;
        double width=(end-begin)/n;

        for(double i=begin;i<=end;i+=width){
            sum+=fx.func(i);
            if(!(i==begin||i==end))sum+=fx.func(i);
        }
        sum*=width/2;
        return sum;
    }

    /**
    * 
    * @param fx    :被積分関数
    * @param begin :積分下限
    * @param end   :積分上限
    * @param n :帯の数　多ければ多いほどよい近似値が得られる
    * @return シンプソンの公式による定積分の近似値
    */
    public static double simpson_rule(Function fx,double begin,double end,int n) {
        if(begin==end)return 0;
        if(fx.isOddFunction()&&begin==-end)return 0;

        double sum=0;
        double width=(end-begin)/n;

        for (int i = 0; i < n+1; i++) {
            double temp;
            temp=fx.func(begin+i*width);
            if(!(i==0||i==n))temp*= i%2==0?2:4;
            sum+=temp;
        }
        sum*=width/3;
        return sum;
    }
}

一応ディレクトリ構成も…。

/javas
   |___lab
   |    |__Main.java
   |
   |___math
         |__Function.java
         |
         |__Integral.java

ちなみにlabっていうパッケージ名には特に意味はないよ。なんとなく実験的なのするみたいなニュアンスはあるかもだけど。

そこまでソースコードで難しい挙動をしているところはないと思うけど、二点だけ補足をするなら、

$${\int_a^a f(x)dx=0}$$

被積分関数が奇関数のとき、

$${\int_{-a}^{a}f(x)dx=0}$$

という性質があるよ。これは、グラフを見てみたらわかるかもしれないね。

なお、いまは$${\int_2^6 \frac{dx}{x}}$$の計算をシンプソンの公式でしているけど、計算法方法や被積分関数、積分範囲も簡単に変えることができるよ。

その場合は、Main.javaのsimpson_rule()をtrapezoidal_rule()に変えたり、func()の中身を変えたり、数値積分の関数の引数の(…2,6…)を変えたりすればいいよ。

また、計算の精度を上げるには最後の引数の値を大きくすればいいよ。

終わりに

ミケ：ふぅ～、今日も疲れた～。

マロ：たっだいま～！

ミケ：おー、マロおかえり！

マロ：どうしたの、ぐったりしてさ。

ミケ：実はさ、数値積分のプログラムを作ったんだ！

マロ：ほうほう、それは面白そうじゃないか。僕にも教えてよ！

ミケ：いいよ～、まずはね～…。

飼い主さんが返ってくるまでの数時間、二人は数学の話で盛り上がったのでした。

---

　　　　　　　誤字脱字、質問ありましたら遠慮なくコメントをください。
　　　　　　　　　　　　　　　　　　　　　　　　　　　　　by ミケ