[Pythonによる科学・技術計算] モンテカルロ積分, 数値計算, numpy

#はじめに
モンテカルロ積分法により
$ I = \int_0^1 \frac{4}{1+x^2} dx = \pi $を計算する。

#内容
(1) 単純モンテカルロ法による計算。生成された乱数が積分領域にいくつ入るか(多数のボールを投げてマトにいくつ入るか)により面積を評価する。

(2) 平均値法によるモンテカルロ積分計算。関数の定義域を一様にサンプルし， 生成されたx点におけるy値の平均値を用いる方法。一般に，yの平均値を$y_{av}$とすると，

積分$I= \int_a^b f(x) dx \sim \frac{b-a}{N}y_{av}$

となる。そこで$y_{av}$をモンテカルロ法で決定しようというのがこの方法のミソである**[1]**。

(3)**平均値法によりる多重積分計算。**例として，
$\int_1^2 \int_1^2 \frac{1}{x+y} dx =10 ln2-6ln3 \sim 0.33979807 $
を計算する。

---

#コード(1): 単純モンテカルロ法による計算

simple_Monte_Carlo.py

import numpy as np
from random import random
"""
単純モンテカルロ積分法: 試行回数 Nを10から10^7まで変化させる
"""
def f(x):
    return 1.0/(1.0+x**2) # 関数の定義

N_calc_list = [10, 10**2, 10**3, 10**4, 10**5, 10**6, 10**7]

for N in  N_calc_list:
    count = 0.0
    for i in range(N):
        x = random()  # [0,1]までの一様乱数をxに格納
        y = random()  # [0,1]までの一様乱数をyに格納
        if y < f (x):   #もし”マト”に入ったらそれをカウントする
            count +=1.0
    area = 4*count/N # 積分結果。 π の評価

    print(N, ", ", area, ", ", abs((np.pi-area)/np.pi)) #結果の出力

---

#結果(1)

10 ,  4.0 ,  0.2732395447351627
100 ,  3.08 ,  0.01960555055392467
1000 ,  3.184 ,  0.01349867760918959
10000 ,  3.12 ,  0.006873155106573032
100000 ,  3.14716 ,  0.0017721414021786754
1000000 ,  3.143488 ,  0.0006033075001118286
10000000 ,  3.1414272 ,  5.266551333578959e-05

左から， モンテカルロ法で用いた乱数の総数N, 積分値, 相対誤差(πからの誤差)

---

#コード(2):平均値法によるモンテカルロ積分計算

The_mean_value_method_Monte_Carlo.py

import numpy as np
from numpy.random import rand
"""
平均値法によるモンテカルロ積分
The mean value method
"""

    
N_calc_list = [10, 10**2, 10**3, 10**4, 10**5, 10**6, 10**7]
for N in  N_calc_list:
    x_array=rand(N)
    y_array=4.0/(1.0+x_array**2)
    area = (1.0-0.0)*np.sum(y_array)/(N) # y_avとして， Σ_i yi/Nを計算し， それに積分区間(1-0=1)を掛けている

    print(N, ", ", area, ", ",  abs((np.pi-area)/np.pi))#結果の出力

---

#結果(2)

10 ,  3.18542467193 ,  0.0139521647705
100 ,  3.03821388912 ,  0.0329064827523
1000 ,  3.14697964989 ,  0.0017147341794
10000 ,  3.14560900784 ,  0.00127844526391
100000 ,  3.14380423975 ,  0.000703969738006
1000000 ,  3.14195518509 ,  0.000115397359081
10000000 ,  3.14140220827 ,  6.06206294417e-05

左から， モンテカルロ法で用いた乱数の総数N, 積分値, 相対誤差(πからの誤差)

---

#コード(3):平均値法による二重積分

double_integral.py

import numpy as np
from numpy.random import rand
"""
平均値法によるモンテカルロ積分(二重積分)
"""


N_calc_list = [10, 10**2, 10**3, 10**4, 10**5, 10**6, 10**7]
for N in  N_calc_list:
    x_array=rand(N)+1.0 # 区間[1,2]の一様乱数列をx_arrayへ格納
    y_array=rand(N)+1.0 # 区間[1,2]の一様乱数列をy_arrayへ格納
    z_array=1.0/(x_array+y_array)  # z=1/(x+y)の列をz_arrayへ格納
    area = (2.0-1.0)*(2.0-1.0)*np.sum(z_array)/(N)  #積分計算
    print(N, ", ", area) #結果の出力

---

#結果(3)

10 ,  0.346923895691
100 ,  0.343570822229
1000 ,  0.339161537421
10000 ,  0.340241114706
100000 ,  0.339920052165
1000000 ,  0.339757274305
10000000 ,  0.33977097358

左から， モンテカルロ法で用いた乱数の総数N, 積分値。厳密解は$10 ln2-6ln3 \sim 0.33979807 $である。

---

#参考文献
[1] Mark Newman, "Computational Physics", Chap10, Createspace Independent Publishing Platform (2012)

---

A comment to non-Japanese persons:

The two codes above described are both simple for performing numerical integrations using the Monte Carlo methods. One can easily understand how the codes work. The detailed information is accessible via ref. [1].