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

最新コメント
リンク

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