プロフィール

kazutabakun

Author:kazutabakun
好奇心旺盛でいろいろなことに興味を持つタイプ。

最新記事
Amazon
時計
最新コメント
最新トラックバック
カテゴリ
月別アーカイブ
FC2ブログランキング
気に入ったら下記をクリックしてください。

FC2Blog Ranking

検索フォーム
RSSリンクの表示
リンク
ブロとも申請フォーム

この人とブロともになる

スポンサーサイト

上記の広告は1ヶ月以上更新のないブログに表示されています。
新しい記事を書く事で広告が消せます。

EXCELで遊ぼ(1): 円周率の計算

この週末に少し時間があったので、ひょっと思いついて円周率の計算でもしてみようかと思いました。
真っ当な計算は既にいろいろ披露されています。この手の計算は、スーパーコンピュータのベンチマーク
にも使われています。また歴史を見るとアルキメデスが内接/外接する多角形から円周率を計算したのが
初期の頃の試行として有名です。

450px-Domenico-Fetti_Archimedes_1620.jpg 
アルキメデス、シシリア島のシラクサに住んでいました。
ここでアルキメデスの原理(比重)やその他革新的な研究(思索)を行っていました。
第2ポエニ戦争(ハンニバルとローマの戦争)でシラクサがローマ軍に攻められたときに、
アルキメデスの考案した新兵器がローマ軍を悩ましたと記録に残っているそうです。
ローマ軍司令官が命は奪うなと厳命したのですが、地面に図を書きながら
一生懸命考えていたアルキメデスが、一兵士の言うことをきかなかったために殺害されたてしまいました。

=======
さてここで使う計算の方法は、確率を使って円周率(pai)を計算するという物です。
もっともこの方法もモンテカルロ法により計算としていろいろな所に掲載されていますので、
あまりオリジナリティがある代物ではないのですが、今回エクセルを使って
簡便に計算してみることにトライしました。

さて、下記の図のように正方形に内接する円を考えます。

円周率図

この座標のXとYがプラスになる領域を考えます。中心を原点として右上のみ。(第一象限)
そのときの内接する円弧で囲まれる面積S1は下記の式でしめされます。

S1=円周率x半径x半径/4

つまり円の面積を4分割した領域になります。その背景の正方形の面積をS2とすると
半径x半径になりますので円周率は下記の式で計算することができます。

円周率= (S1 x 4) / S2

ここで次の仮定を導入します。
1) ものすごく沢山矢を放つと半径を1辺とした正方形内にその痕が一様に分布する。
2) 矢を射るとその合計がS2本とした場合そのうちこの円弧に入る矢の数がS1にほぼ等しくなる。
これは上記(1)の仮定から円弧内に入る矢がS1/S2の確率で決まるというのが根拠です。
3) したがってS1本、S2本にした数を上記の式に入れると円周率が計算できる。

はたしてどれぐらいの精度でこのことが言えるのかエクセルを使って実際に計算してみました。

一辺を1にした場合を考えました。理由はエクセルでは0から1までの乱数を発生させる
=RAND()という関数を使うためです。また精度を出すためには、相当たくさんの試行が必要です。

エクセルでは、X=RAND(), Y=RAND()にするとX軸、Y軸に点P(X,Y)が0<X<1、0<Y<1の範囲で
ランダムに生成されます。この数が多いほど領域が塗りつぶされます。

内接する円弧をあらわす式はX^2+Y^2=1になりますのでこの足し算の結果が1より小さいとその点Pは
円弧の中に入ることになります。なおここで記号X^2はXの2乗を意味します。
さて円弧の中に入る点の特徴ですが、下図を参照してください。
ここでPinは円弧内の点をあらわしていますが、中心からの距離はあきらかにRよりは
短いことが分かります。
すなわち円弧内の点のR(R=1とした場合)については、
xの2乗とYの2乗の和は1より小さいということになります。
円弧の内部点

総生成数をS2個とした場合、円弧の中に入るPの合計数をS1個と考えます。
上記のような考え方を使って計算しました。その計算結果の一部を下記表で示します。
ここでRound Pというのは円弧に入った場合は1、外れた場合は0になるようにしています。
Square Pは全て1であとで合計生成数に使うためです。

円周率の計算

一つのシートで10,000回の試行ができるようにしました。あまり多い回数にすると
エクセルの動きが遅くなりますので。
その試行を30回程度繰り返し合計本数を300,000個にした場合のS1の数を求めて円周率を計算しました。
円周率=3.14...程度までの結果が出ましたが、小数点3桁以下はまだ誤差がありました。その原因は
1)試行回数がまだ少ない
2)エクセルで発生させるランダム関数が本当にランダムか?
等があります。

円周率の計算結果

1,000回乱数を発生させた円弧の中に入る点をプロットして図示できるようにしました。
何となく円弧の形になっていますが、スキマが多くこのぐらいの回数だと完全に面積と
同じになるとは言えない状態です。

円周率確率分布

スポンサーサイト

Comment

非公開コメント

上記広告は1ヶ月以上更新のないブログに表示されています。新しい記事を書くことで広告を消せます。