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

最新コメント
リンク

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