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