snovaのブログ

プログラミングとか、日常のこととか、アウトプットしたほうがよいと聞いたので

Pythonでざっくりと数値解析の基礎をしてみる(3_数値積分編)



【スポンサーリンク】



イントロダクション

今回は数値積分を解きます。
区分求積法、台形公式、モンテカルロ積分の紹介をします。

目次

区分求積法

長方形を並べていき、その面積の合計を積分値として出す方法です。
非常に単純ですが、高校数学では積分の概念を学ぶときに出てくるほど重要な方法です。

幅を小さくしていくと

のように、誤差が少なくなっていきます。

では、実際に $$ y = \int_{0}^{\pi} \sin x dx $$ を解いてみます。

ちなみに、答えは2です。

以下、プログラム

import numpy as np
import matplotlib.pyplot as plt

def sectional_measurement():
    w = 1 # 長方形の幅
    x = np.arange(0, np.pi, w) # x
    y = np.sin(x) # y

    val = np.sum(y * w) # 長方形の面積の合計
    print(val)

if __name__ == '__main__':
    sectional_measurement()

結果

1.8918884196934453

長方形の幅を w = 0.1 のように小さくすると

1.9995479597125976

のように、誤差が小さくなります。

台形公式

長方形の代わりに台形を使います。

簡単にグラフで表すと

刻みの幅を小さくすると

長方形より誤差が少なくなります。

同様に、 $$ y = \int_{0}^{\pi} \sin x dx $$ を解いてみます。

import numpy as np
import matplotlib.pyplot as plt

def trapezoid():
    w = 1 # 台形の高さ(刻み幅)
    x = np.arange(0, np.pi, w) # x
    y1 = np.sin(x) # 台形の上底
    y2 = np.roll(y1, 1) # 台形の下底
    y3 = y1 + y2 # 上底+下底
    y3 = y3[1:] # rollを使ったので、0番目を除く
    y = y3 * w / 2 # (上底+下底)*高さ/2

    val = np.sum(y) # 台形の合計
    print(val)

if __name__ == '__main__':
    trapezoid()

結果 :

1.8213284156635117

台形の高さを w = 0.1 のように小さくすると

1.997468926590933

区分求積法より誤差が大きくなったのはなぜでしょう?

ちなみに、 scipy には台形公式のモジュールがあります。

import numpy as np
import matplotlib.pyplot as plt
from scipy import integrate

def trapezoid():
    w = 0.1 # 台形の高さ(刻み幅)
    x = np.arange(0, np.pi, w)
    y = np.sin(x)
    int_scipy = integrate.cumtrapz(y, x) # 積分の計算
    print(int_scipy[-1])

if __name__ == '__main__':
    trapezoid()

結果 :

1.9974689265909336

さらに進化させて、2次式で近似するとシンプソン公式になります。

モンテカルロ積分

期待値を使った数値積分です。
プログラムでは、乱数によって生成された値を使い計算された値が、積分される領域内にあるかを判定し、数えることで、積分を行います。
メリットは複雑な積分計算でもできるところです。

サンプル数を増やすと精度が上がります。

同様に、 $$ y = \int_{0}^{\pi} \sin x dx $$ を解いてみます。

import numpy as np
import matplotlib.pyplot as plt

def Monte_Carlo():
    N = 10000
    np.random.seed(0)
    x = np.pi * np.random.rand(N) # xを[0, pi]の範囲からランダムにN個選択
    y = 2 * np.random.rand(N) - 1 # yを[-1, 1]の範囲からランダムにN個選択
    f = np.sin(x) # もとの関数

    count_y = np.sum((0 <= y) & (y <= f)) # yが[0, f(x)]の範囲にあるか判定し個数を計算
    val = 2*np.pi * (count_y / N) # 積分値を計算
    print(val)

if __name__ == '__main__':
    Monte_Carlo()

結果 :

1.9603538158400309

半径2の円の体積も出してみました。
答えは $$ V = \frac{4}{3} \pi r^3 = 33.5103216... $$

import numpy as np
import matplotlib.pyplot as plt

def Monte_Carlo_sphere():
    N = 10**5
    np.random.seed(0)

    r = 2 # 球の半径
    x_start = -r # xの範囲の最小値
    x_stop = r # xの範囲の最大値
    x = (x_stop - x_start) * np.random.rand(N) + x_start

    y_start = -r # yの範囲の最小値
    y_stop = r # yの範囲の最大値
    y = (y_stop - y_start) * np.random.rand(N) + y_start

    z_start = -r # zの範囲の最小値
    z_stop = r # zの範囲の最大値
    z = (z_stop - z_start) * np.random.rand(N) + z_start

    f = np.sqrt(x**2 + y**2 + z**2)
    count_f = np.sum(f <= r) # 球内部に乱数で取った値があるのかの判定
    int_val = (2*r)**3 * (count_f / N) # 積分値の計算
    print(int_val)
    print(4*np.pi*(r**3) / 3) # ついでに公式から導いた式も出しておく

if __name__ == '__main__':
    Monte_Carlo_sphere()

結果 :

33.39072
33.510321638291124

まとめ

数値積分を台形公式、モンテカルロ積分で実践しました。
どのくらいの精度が欲しいかにもよりますが、高速に計算できます。

参考文献