今回は「モンテカルロ法」で円周率を求めてみたいとおもいます。
モンテカルロ法は乱数を使ったシミュレーション方法です。
そのため今回はfortranのサブルーチンの一つである、random_numberを使います。
ではさっそくやっていきましょう。
モンテカルロ法で円周率を求める
円周率は以下のような方法で考えます。
まず、一辺2の正方形の中に半径1の円(円の面積は\(\pi\))を考えます。
次にその中にランダムに点を打っていきます(\(N\)個)。
そして円の中にあたった点だけを数えます(\(N_h\)個)。
すると円の面積と正方形の面積の比は
$$
円の面積 : N_h = 正方形の面積 : N
$$
つまり、
$$
円の面積(=\pi) = 正方形の面積 \times \frac{N_h}{N}\tag{1}
$$
となり、円周率が求まります。
これをfortranで書いていきます。
fortranで数値計算
fortranで書くと以下のようになります。
program montecarlo
implicit none
integer :: N,Nh,i
real :: x,y,S
real,parameter :: pi=acos(-1.0)
Nh = 0
N = 80000
do i = 1,N
call random_number(x)
call random_number(y)
x = (1 - (-1))*x -1
y = (1 - (-1))*y -1
if(x**2 + y**2 < 1)then
Nh = Nh +1
endif
enddo
S = 4.0 * Nh/N
write(*,*) "S=",S
write(*,*) "pi=",pi
end
以下で解説します。
\(x,y\)を-1~1の乱数として作成
サブルーチンを使い\(x,y\)を乱数にします。
さらに一辺\(2\)正方形の中に点を打ちたいので\(x,y\)を\(-1 < x,y < 1\)に限定します。
\(a<x<b\)の範囲に限定したいときは
$$
x = (b – a) \times x + a
$$
とすればよいので
call random_number(x) call random_number(y) x = (1 - (-1))*x -1 y = (1 - (-1))*y -1
とします。
これで正方形内にランダムに点を打つことができました。
円の中にあたった点の数を数える
if文を使い、円の内部(\(x^2 + y^2 < 1\))に\(x,y\)がある数をカウントします。
つまり、この条件を満たすたび、\(N_h\)を増やします。
すなわち
if(x**2 + y**2 < 1)then Nh = Nh +1 endif
となります。
面積計算
上の二つのステップを繰返しループさせ、できるだけ多くの点を打ちます(今回は\(N=80000回)\)。
そして(1)式を使って円の面積(=\(\pi\))を計算します。
S = 4.0 * Nh/N
すると結果は
S= 3.14000010 pi= 3.14159274
くらいになり、だいたい近い値が出ました！！
このようにしてだいたいの円周率は求めることができますが、工夫次第でより精度の高い値を出すこともできます。
興味があればぜひやってみてください。
