Scilabで円周率の近似値

3月14日は円周率の日だからというわけではないですが、Scilabで楽しむ確率論に倣ってモンテカルロ法を用いた円周率の近似値計算を行いました。

001_20131105060858590.jpg

Fig.1: πのパイ(Wikipedia Pi pie2.jpg)


半径1の円の面積はπなので、その四分の一の面積はπ/4です。
したがって0≦x≦1, 0≦y≦1の範囲にランダムに点を打ち、半径r=1の円の中に入っている数を数えれば、円周率の近似値を求めることが出来ます。

今回はn=10000個の乱数の組(x,y)が半径r=1の円の中に入る数から円周率を求め、更にこれをm=1000回繰り返したときどの程度のばらつきがあるのかをヒストグラムにしました。


モンテカルロ法による円周率計算


0≦x≦1, 0≦y≦1の領域に一様分布に従うように点を打ったとき、半径r=1の円の中に入る点の数がπ/4に比例します。
以下のようなScilabプログラムにて円周率の近似値を得ることが出来ます。

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


002_20131105060857b76.png

Fig.2: [0,1]×[0,1]の領域内にランダムに点を打ったとき、半径1の領域に打たれる点の頻度が近似的にπ/4になる。


繰り返し計算


以上の計算を更にm=1000回繰り返してどの程度のばらつきがあるのかをヒストグラムにしました。
乱数を格納するために大きなメモリを確保する必要があります。これにはstacksize関数を利用します。

m=1のときの計算も含めたScilabプログラムはcircle_sce.txtとなりました。

003_2013110506085769a.png

Fig.3: (n,m)=(10000,1000)のときのヒストグラム


関連エントリ




参考URL




付録


このエントリで使用したScilabのシミュレーション用ファイルを添付します。ファイル名末尾の".txt"を削除して、"_"を"."に変更すれば使えるはずです。(参考:ねがてぃぶろぐの付録)


参考文献/使用機器




フィードバック



にほんブログ村 その他趣味ブログ 電子工作へ

 ↑ 電子工作ブログランキング参加中です。1クリックお願いします。


コメント・トラックバックも歓迎です。 ↓      


 ↓ この記事が面白かった方は「拍手」をお願いします。

tag: Scilab モンテカルロ解析 乱数 ヒストグラム 円周率 

FC2カウンター
カテゴリ
ユーザータグ

LTspiceAkaiKKRmachikaneyamaScilabKKRPSoCOPアンプCPAPIC強磁性モンテカルロ解析常微分方程式トランジスタode状態密度インターフェースecalj定電流スイッチング回路PDS5022DOS半導体乱数シェルスクリプトレベルシフトHP6632Aブレッドボード分散関係温度解析トランジスタ技術R6452A可変抵抗I2Cセミナー確率論反強磁性非線形方程式ソルバ絶縁偏微分方程式バンド構造熱設計数値積分バンドギャップカオスA/DコンバータフォトカプラシュミットトリガGW近似LEDLM358ISO-I2C三端子レギュレータ数値微分サーボ直流動作点解析カレントミラーマフィンティン半径TL431PC817C発振回路74HC4053USBアナログスイッチbzqltyFFTチョッパアンプ2ちゃんねる補間量子力学開発環境電子負荷標準ロジックパラメトリック解析アセンブラ基本並進ベクトルブラべ格子単振り子BSchLDAイジング模型繰り返しMaximaキュリー温度位相図状態方程式失敗談スピン軌道相互作用六方最密充填構造相対論FET抵抗コバルト不規則合金TLP621ewidthGGAQSGWgfortranランダムウォークラプラス方程式スイッチト・キャパシタcygwin熱伝導SMPスレーターポーリング曲線三角波格子比熱LM555条件分岐TLP552MCUNE555UPSTLP521QNAPマントルテスタFXA-7020ZR過渡解析詰め回路ガイガー管ダイヤモンド自動計測Writer509データロガー固有値問題VESTAスーパーセルOpenMP差し込みグラフ平均場近似起電力awk仮想結晶近似VCAubuntufsolveブラウン運動熱力学第一原理計算井戸型ポテンシャルシュレディンガー方程式面心立方構造fccウィグナーザイツ胞interp12SC1815L10構造非線型方程式ソルバFSMキーボードTeX結晶磁気異方性初期値OPA2277化学反応等高線ジバニャン方程式ヒストグラム確率論三次元フィルタRealforcePGAフェルミ面正規分布固定スピンモーメント全エネルギースワップ領域リジッドバンド模型edeltquantumESPRESSOルチル構造岩塩構造二相共存ZnOウルツ鉱構造BaOフォノンデバイ模型multiplotgnuplotc/aノコギリ波合金クーロン散乱ハーフメタル半金属CapSenseマンデルブロ集合マテリアルデザインSICGimpCK1026MAS830L円周率トランスPIC16F785凡例線種シンボルLMC662ヒストグラム不規則局所モーメント文字列疎行列TS-110TS-112Excel直流解析等価回路モデル入出力トラックボールPC軸ラベルAACircuitP-10フラクタル境界条件連立一次方程式Ubuntuifortパラメータ・モデルspecx.f関数フィッティング最小二乗法Crank-Nicolson法陰解法日本語EAGLEMBEグラフの分割負帰還安定性ナイキスト線図熱拡散方程式HiLAPW両対数グラフ片対数グラフ縮退

最新コメント
リンク

にほんブログ村 その他趣味ブログ 電子工作へ