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カウンター
カテゴリ
ユーザータグ

LTspiceAkaiKKRmachikaneyamaScilabKKRPSoC強磁性CPAPICOPアンプecalj常微分方程式モンテカルロ解析状態密度トランジスタodeDOSインターフェース定電流スイッチング回路PDS5022半導体シェルスクリプト分散関係レベルシフト乱数HP6632AR6452A可変抵抗トランジスタ技術温度解析ブレッドボードI2C反強磁性確率論数値積分セミナーバンドギャップバンド構造偏微分方程式非線形方程式ソルバ熱設計絶縁ISO-I2Cカオス三端子レギュレータLM358GW近似マフィンティン半径A/DコンバータフォトカプラシュミットトリガLEDPC817C発振回路数値微分直流動作点解析サーボカレントミラーTL431アナログスイッチUSB74HC4053bzqltyVESTA補間電子負荷アセンブライジング模型BSch量子力学単振り子2ちゃんねるチョッパアンプLDA開発環境基本並進ベクトルFFT標準ロジックブラべ格子パラメトリック解析抵抗SMPMaxima失敗談ラプラス方程式繰り返し位相図スイッチト・キャパシタ熱伝導状態方程式キュリー温度gfortranコバルトTLP621不規則合金Quantum_ESPRESSO六方最密充填構造ランダムウォーク相対論ewidthスピン軌道相互作用FETQSGWVCAcygwinスレーターポーリング曲線GGA仮想結晶近似PWscfシュレディンガー方程式LM555ハーフメタル固有値問題NE555最小値ガイガー管QNAPUPS自動計測ダイヤモンドマントルTLP552格子比熱最適化MCU井戸型ポテンシャル最大値xcrysdenCIF条件分岐詰め回路フェルミ面差し込みグラフスーパーセルfsolveブラウン運動awk過渡解析起電力三角波第一原理計算FXA-7020ZRWriter509Ubuntuテスタ熱力学データロガーTLP521OpenMPubuntu平均場近似MAS830LトランスCK1026PIC16F785PGA2SC1815EAGLEノコギリ波負帰還安定性ナイキスト線図MBEOPA2277P-10フィルタCapSenseAACircuitLMC662文字列固定スピンモーメントFSMTeX結晶磁気異方性全エネルギーc/a合金multiplotgnuplot非線型方程式ソルバL10構造正規分布等高線ジバニャン方程式初期値interp1fcc面心立方構造ウィグナーザイツ胞半金属デバイ模型電荷密度重積分SIC二相共存磁気モーメント不純物問題PWgui擬ポテンシャルゼーベック係数ZnOウルツ鉱構造edeltquantumESPRESSOフォノンリジッドバンド模型スワップ領域BaO岩塩構造ルチル構造ヒストグラム確率論マテリアルデザインフラクタルマンデルブロ集合キーボードRealforceクーロン散乱三次元疎行列縮退化学反応関数フィッティング最小二乗法Excel直流解析PCTS-110TS-112日本語パラメータ・モデル等価回路モデルcif2cell入出力陰解法熱拡散方程式HiLAPW両対数グラフCrank-Nicolson法連立一次方程式specx.fifort境界条件片対数グラフグラフの分割円周率ヒストグラム不規則局所モーメントGimpシンボル軸ラベル凡例線種トラックボール

最新コメント
リンク

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