【モンテカルロ法による円周率の計算 】
MATLABで次のように入力してみましょう.
n = 10;
x = rand(1,n)
y = rand(1,n)
r = x.^2 + y.^2
p = sum(r <= 1)/n*4
plot(x,y,'x');
hold on;
ezplot('cos(t)','sin(t)',[0 pi/2]);
hold off;
何度か繰り返すと異なる値が得られます.
Scilabではグラフのための命令が異なります.
plot2d(x',y',-2);
ezplot('cos(t)','sin(t)',[0 pi/2]);
nを変えて何度も計算するために
次のような関数を用意します.ファイル名は
"monte.m"とします.
function p = monte(n)
x = rand(1,n);
y = rand(1,n);
r = x.^2 + y.^2;
p = sum(r <= 1)/n*4;
MATLAB上で
>> monte(100)
のようにします.100を1000や10000に変えると結果は少しずつ
円周率の値に近づいていきます.引数に同じ値を与えても
乱数が異なるため実行するたびに結果が異なります.
同じ乱数を生成するには
rand('seed',0)
のようにします.0を他の整数に変えると乱数のパターンが変わります.