clear; // *** 計算のための設定 *** stacksize('max') // スタックサイズを最大値に設定 // ********************************************************* // // n個の乱数(x,y)の組を使った1回(m = 1)の計算 // // ********************************************************* n = 10000; x = rand(1,n); y = rand(1,n); // グラフの縦横比を等しくする h = scf(0); // ウィンドウを作成 clf(); // 表示をクリア ha = h.children(1); // Axes(座標軸)オブジェクトへのハンドルを取得 ha.isoview = "on"; // 座標軸の縦横比を等しくする ha.data_bounds = [0, 0; 1, 1]; // 座標軸表示範囲の設定 // 乱数の組(x,y)をプロット plot(x,y,'.r','markersize',1); // しきい値となる半径1の円のプロット xx = linspace(0,1); plot(xx,sqrt(1 - xx .^ 2),'-b'); // グラフのラベル xlabel("x"); ylabel("y"); // 1回の計算から得られた円周率 p1 = sum(bool2s(x ^ 2 + y ^ 2 <= 1.0)) / n * 4 // ********************************************************* // // n個の乱数の組(x,y)を使ったm回の計算の平均とヒストグラム // // ********************************************************* m = 1000; X = rand(m,n); Y = rand(m,n); P = sum(bool2s(X .^ 2 + Y .^ 2 <= 1.0),"c") ./ n .* 4; scf(1); clf(); histplot(100,P); // m回の計算から得られた円周率 pm = mean(P) // 円周率の文献値 %pi