algorithm第4回

file://H:\_サーバー\lachesis20111116\Manual\Algorithm\algorithm4.html
アルゴリズム(数値積分編)
[ モンテカルロ法( Monte Carlo Method ) ] [ 課題 ]
数値積分
3. モンテカルロ法( Monte Carlo Method )
各数の前後に何の因果関係も存在しない数列を乱数 ( random number ) という。
ここでは、各数が独立で等しい確率で現れる一様分布 ( uniform distribution ) 乱数を用いる。
すなわち、各数の前後に何ら因果関係も存在しないが、各数ごとに出現確率をとれば、ほとんど同じ数だけ出現する。
モンテカルロ法( Monte Carlo Method )は確率分布乱数を用いたシミュレーションの総称であり、確率的な因子を含む問題はすべてモンテカルロ法の対象である。
ここでは、
をモンテカルロ法を用いて近似的に求める方法を紹介する。この積分値は図 3.4 で示される斜線部分の面積である。
図 3.4
の意味
図 3.5 で示されるように、斜線部分を完全に含む最も小さい四角形を考え、その面積 V を考える。次に一様分布乱数を用いてこの四角形内に点をランダムにばらまく。
図 3.5 モンテカルロ法による
の近似解
つまり、点 vt の x 座標 xt と y 座標 ytをそれぞれ
xt = ( b - a ) ・ RND( 1 ) + a
yt = c・RND( 1 )
として点を設ける。ただし、RND( 1 )は区間( 0, 1 )の一様分布乱数であり、BASIC言語の組み込み関数である。四角形内にばらまかれた点のうち、斜線部分に入った点の数 k を数え、乱数の性質 V : S = n : k を用いて
S を求める。
ここで、点vt が斜線部分に入るための必要十分条件は
f ( xt ) ≧ yt
である。これを図 3.6に示す。
1/3 ページ
file://H:\_サーバー\lachesis20111116\Manual\Algorithm\algorithm4.html
図 3.6 点vt が斜線部分に入るための必要十分条件
図 3.7 と図 3.8に
の近似解を求めるフローチャートとプログラムを示す。また、図 3.8 の実行例は
図 3.7
10 '* Monte Carlo method
20 CLS
30 INPUT A, B, C
40 INPUT N
50 K=0
60 V=C*(B-A)
70 FOR I=1 TO N
2/3 ページ
の近似解であり、Nの値を大きくするほど理論値 π/4 に近づくことがわかる。
の近似解を求めるフローチャート
file://H:\_サーバー\lachesis20111116\Manual\Algorithm\algorithm4.html
80 X=(B-A)*RND(1)+A
90 Y=C*RND(1)
100 IF 1/(1+X^2) >= Y THEN K=K+1
110 NEXT I
120 S=V*K/N
130 PRINT "* A=";A, "B =";B, "C =";C
140 PRINT "* 全体にばらまいた点の数 N =";N
150 PRINT "* 面積は";S, "4倍すると、π=";4*S
160 END
図 3.8
の近似解を求めるプログラム
[ 課題 ]
1. 変数 x の変域 [a, b] における関数 y=f (x) を表示するプログラムを作成せよ。
ヒント
i. x の変域内を500等分し、x0, x1, …, xnの各 xi に対する yi を求め、配列に格納する。
ii. 変域内における関数の最大値、最小値を求める。つまり、先の配列中から最大値 MAX と最小値 MIN を求める。
iii. 一般に、直接画面に点(xi, yi) を表示することはできないので(BASICにおいてy軸方向に打てる点の数は一般に400点である)、yiの値が画面に適度に収まるように調整する。図 3.9参照。
y軸方向に打てる点の数を300とすると、MIN から MAX まで変化する yi の値を変換する必要がある。
yi=300*( yi - MIN ) /(MAX-MIN)
iv. そして、求まった点 ( xi, yi )を画面上にプロットするのだが、x軸を500等分、y軸を300等分したのは原点を ( 100, 350 )にしてもらうためです。一般のグラフのy軸の方向とコンピュータ画面のy軸方向が異なっている
ことに気をつけて、求めた各点をプロットして下さい。
図 3.8
2. 楕円を表示するプログラムを作成せよ。
楕円上の各点は( x, y )は、( a * cos(θ), b * sin(θ))とし、θを0∼360まで変化させよ。(プログラミングにおいてはθという記号は使わないで"th"、とでもしておいて下さい。)
また、係数 a, bもθの関数とすることで面白い模様を表示することができる。
例、a=b=40 * sin( 8 * θ) + 50
3/3 ページ