突然ですが、初めて円周率を自分で計算したときのことを覚えていますか?私は鮮明に覚えています。そう、あれば小学校5年の算数の授業でした。
先生に急にダンボールと糸を渡され、「円周率3.14を導け」と言われたのです。当時はうまく理解できていなかった気がしますがなんとか手伝ってもらいながら無事3.14に近い値を算出したことを覚えています。そのときの手順は次の通り。ものすごくアナログです(笑)
1. コンパスでダンボールに円を描く
2. 円を切り抜く
3. 円の周りに沿って糸を一周巻いて切る
4. 糸の長さを測る
5. 糸の長さをダンボールの円の直径で割る
ちなみに手順3が一番難しいです。適当に切ると長すぎたり短くなったりするんです。コツはちょっと余分に巻いて重なっている部分をカッターで一気に切るんです!こうするとぴったり切れるんです。
もし暇すぎで死んでしまいそうだったらやってみてください。円周は直径に円周率をかけることで求められるので、円周率をもとめるには円周(糸の長さ)をダンボールの円の直径で割れば円周率っぽいものが得られるというものです。
\(\large r \): 半径
\(\large \pi\): 円周率
\(\large l \):円周
上のように設定したとき、円の円周は次の式で求めることが可能です
\(\large 2\pi r = l\)
これを\(\large\pi\)で整理します
\(\large\pi = \frac{l}{2r}\)
大数の法則で近似値を求める方法
私もあれから大きくなりました。いまではもうおっさんです。おっさんがダンボールと糸でちまちまやっててもかっこよくないのでプログラムで円周率を求めてみましょう。
その前に問題です。
1辺が1の正方形があります。
この四角の半分に色を塗ります。
この四角形の上にランダムに点を打つと青い部分に入る確率は何%でしょうか?正解は、50% です。
もちろん次のように50%にならない場合もありますが点を増やすに従って結果は\(\large\frac{1}{2}\)に近づいていきます。いわゆる大数の法則というものです。
では、なぜ\(\large\frac{1}{2}\)に近づくか。答えは感覚的にも分かる通り、青い部分の面積が全体のちょうど50%だからです。このことからわかるように面積の比率と穴が打たれる確率は一致するものと考えることが可能です。この考え方で円周率を求めていきます。
円周率を近似する方法
円周率の近似は次の図形を元に考えます。
この上にランダムな点をたくさん打ち込みます。
わかりやすいように円の内側に打たれた点と外側に打たれた点を
区別しておきましょう。
正方形に含まれる扇型の円の面積は半径1の円の\(\large\frac{1}{4}\)の大きさなので次のようになります。
\(\large \frac{1\times1\times\pi}{4} = \frac{\pi}{4}\)
つまり点をたくさん打ったとき、扇型の内側に点が入る確率は徐々に\(\large\frac{\pi}{4}\)になっていくということです。式にしてみます。全部の点を\(N\)、扇型に入った点を\(P\)とすると次のように書くことができます。
\(\large\frac{P}{N} = \frac{\pi}{4}\)
これを\(\large \pi\)で整理すると
\(\large\pi = \frac{4P}{N}\)
となります。
\(\large \pi\)を求めるには、扇型に入った点の数\((P\))を全ての点(\(N\))で割って4倍すれば求めることができます。
ということで、この内容をプログラムで実装していきます。ようやくプログラムの話ができますね。
プログラムでπを求めるのに必要なこと。
1辺が1の正方形に無数に点を打つにはどうすれば良いか?0から1までのランダムな数値を2つ用意した
それぞれを\(\large x\)座標、\(\large y\)座標を表すことで正方形の中に収まる点を打つことができます。
扇型の中に打たれた点を判断するにはどうすれば良いか?座標\(\large x\)と\(\large y\)がわかれば、三平方の定理で斜辺の長さが求められます。その長さが円の半径1以下であれば扇型の中に打たれた点と判断することができます。
上のことを踏まえて作成したのが次のプログラムです。
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 |
import random import time point = 0 N = 100 repeat = 5 start = time.time() for count in range(repeat): N = N * 10 point = 0 for i in range(N): x = random.random() y = random.random() if x * x + y * y < 1.0: point += 1 # 円周率の近似解 pi = 4.0 * point / N process_time = time.time() - start print('N=' + str(N) + " \t" + 'pi=' + str(pi) + " \t" + 'time=' + str(process_time)) |
回数によって円周率の精度と処理時間を確認したいため、1000, 10000, 100000, 1000000, 10000000と点の数を変えて円周率を求める計算をしています。
実行結果はこちら
1 2 3 4 5 |
N=1000 pi=3.188 time=0.0013129711151123047 N=10000 pi=3.1612 time=0.015565156936645508 N=100000 pi=3.14284 time=0.14728093147277832 N=1000000 pi=3.138936 time=1.3955268859863281 N=10000000 pi=3.14151 time=13.65818190574646 |
無事円周率っぽい値が出ました。この円の近似方法ははモンテカルロ法といいます。円周率を求める方法は他にもライプニッツの公式ガウス=ルジャンドルのアルゴリズム等々様々なものがありますので、是非プログラムに落とし込んで実行してみてください。それでは、また。