【fortran】モンテカルロ法で円周率求めてみた!

プログラミング

今回は「モンテカルロ法」で円周率を求めてみたいとおもいます。

モンテカルロ法は乱数を使ったシミュレーション方法です。

そのため今回は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 &lt; 1)then
   Nh = Nh +1
endif

となります。

面積計算

上の二つのステップを繰返しループさせ、できるだけ多くの点を打ちます(今回は\(N=80000回)\)。

そして(1)式を使って円の面積(=\(\pi\))を計算します。

S = 4.0 * Nh/N

すると結果は

S=   3.14000010
pi=   3.14159274

くらいになり、だいたい近い値が出ました!!

このようにしてだいたいの円周率は求めることができますが、工夫次第でより精度の高い値を出すこともできます。

興味があればぜひやってみてください。

タイトルとURLをコピーしました