こんにちは！ほけきよです。
皆さん、πを知っていますか？？あの3.14以降無限に続く円周率です。
昔、どこかのお偉いさんが「3.14って中途半端じゃね？ｗｗｗ3にしようぜ」
とかいって一時期円周率が3になりかけました。でもそれは円じゃなくて六角形だからだめです。全然ダメ。
それを受けて「あほか、円周率をちゃんと教えろ」 と主張したのが東大のこの問題*1
めっちゃ単純な問題。でも、東大受験生でさえ 「普段強制的に覚えさせられたπというやつ、どうやったら求められるの???」 と悩んだことでしょう。
また、普段生活してると
「π求めてぇ」
と悩むこともあるでしょう。今日はそんなみなさんに、様々なπの求め方をお教えします。これで、あらゆる状況で求められるようになりますよ！
東大の問題へのアプローチ2つ
もちろん、πの厳密な値を求めることはできません。今でもπの値は日々計算され続けています。 じゃあ、πより少し小さい値で、うまくπの値を近似できる方法を考えよう。 というアプローチです。
多角形で近似
おそらく一番多かったであろう回答が、この多角形近似です
同じ半径であれば、正多角形はすべて円の中に収まります。正方形も正六角形も正八角形も。 なので、それを利用してやりましょう。正六角形は周と直径の比が3であることは簡単にわかるので
- 正六角形よりも多角形
- sinやcosの値が出せそう
な正八角形(もしくは正十二角形)を選びます。
解法はこんな感じです。
tanの逆関数を使う
この問題に関しては、こんな解法もできます！ 高３のときに習いますね！
- 置換積分を使うと、答えにπが現れる
- かつ、上に凸な関数
- かつ、値を代入した時に計算がしやすい
と言えば、そう、
ですね！！
は、ルートがある分、ちと使いにくいのです。
解法は↓のような感じ
無限級数を覚えておく
フーリエ級数を用いる
世の中にはこんな不思議な式があります
これを理解するためには，Fourier級数を知る必要があります。理系の方なら大学1-2年くらいで学びますね。
打ち切り項数との関係はこんな感じ。
N:1 Value:2.4494897 N:10 Value:3.0493616 N:100 Value:3.1320765 N:1000 Value:3.1406381 N:10000 Value:3.1414972 N:100000 Value:3.1415831
フーリエ級数がわかれば、上の式以外にも、例えばこんな式も作れるようになります
分数なら簡単に計算できるし，πも簡単に求められそうですね^^
ラマヌジャン式を使う
無性にπが求めたくなった時も，この無限級数を知っているだけでOK! あの天才ラマヌジャンが導出した式です
美しい式ですね(白目)
めちゃくちゃ収束が早いことが知られているので，n=0,1,2とかをぶち込んでやるだけでそれなりの精度が出るのがいいところ
n = 0, 1での代入結果がこちら
n:0 Value:3.14158504007123751123 n:1 Value:3.14159265359762196468
n=0で、もう良さげ。すごい精度。
- ちょっと複雑で覚えにくい
- 分子分母の値がでっかくなりすぎて計算がそもそも厳しい
のがたまに傷かな？？
コンピュータを使う
モンテカルロサンプリングする
あなたの眼の前にそこそこいいパソコンがあるなら，モンテカルロサンプリングでπを求めましょう！
を，両方とも-1から1の範囲でランダムに選択
を計算
なら +1,
なら何もしない
- N回繰り返して点をばらまく
だった点の数をNで割る
最終的にこの結果を4倍すればPiが求められます
いいところは，回数をこなせばこなすほど精度が上がるところと、事前に初期値設定が必要ないところ。
点を打つほど円がわかりやすくなってくる
悪いところはPCを痛めつけることになること。精度の収束も悪く、計算に時間がかなりかかります。
N:10 Value:3.200000 Time:0.00007 N:100 Value:3.200000 Time:0.00013 N:1000 Value:3.064000 Time:0.00129 N:10000 Value:3.128000 Time:0.01023 N:100000 Value:3.147480 Time:0.09697 N:1000000 Value:3.143044 Time:0.93795 N:10000000 Value:3.141228 Time:8.62200 N:100000000 Value:3.141667 Time:94.17872
無限に時間と計算資源がある人は，試してみましょう！
ガウス=ルジャンドルのアルゴリズムを使う
もっと精度よく効率的に求めたい！！というアナタ！ ガウス=ルジャンドルのアルゴリズムを使いましょう
ガウス＝ルジャンドルのアルゴリズムは円周率を計算する際に用いられる数学の反復計算アルゴリズムである。円周率を計算するものの中では非常に収束が速く、2009年にこの式を用いて2,576,980,370,000桁（約2兆6000億桁）の計算がされた(Wikipediaより)
なんかすごそう…よっぽど複雑なのかと思いきや、アルゴリズムは超簡単(Wikipediaより)
実際にコードを書いてみて動かした結果がこちら
import numpy as np # update algorithm def update(a, b, t, p): new_a = (a+b)/2.0 new_b = np.sqrt(a*b) new_t = t-p*(a-new_a)**2 new_p = 2*p return new_a,new_b,new_t,new_p # initialize a = 1.0 b = 1/np.sqrt(2) t = 0.25 p = 1.0 print("0 : {0:.10f}".format((a+b)**2/(4*t))) # run for i in range(5): a,b,t,p = update(a,b,t,p) print("{0} : {1:.15f}".format(i+1,(a+b)**2/(4*t)))
結果が
0 : 2.9142135624 1 : 3.140579250522169 2 : 3.141592646213543 3 : 3.141592653589794 4 : 3.141592653589794 5 : 3.141592653589794
2回の更新でモンテカルロサンプリングを超えていることがわかります。しかも更新も一瞬！ かなり優秀なアルゴリズムのようです。
実験で求める
ビュフォンの針
もしあなたが針やつまようじを大量に持っているならば、こんな実験をしてみましょう
- 針の長さを測る(
)
- その半分の長さの間隔で平行線を引く(
)
- 引いた平行線の上に針をぶちまける
- 平行線に交わる針の数を数える
これはビュフォンの針問題と言って、針の数をめちゃくちゃ増やすと
となります。
こうするだけで、なんとが求まります。ね、簡単でしょ？？？
単振動
円周率が求めたいときに、バネを見つけたとします。 それはラッキーですね。早速バネの振動する周期を求めましょう！！
図のように、周期にが含まれているので、ばねの振動する時間を求めるだけで、簡単に
が求まります。
注意点は
- 摩擦があると厳密に周期が求められない
- 空気抵抗があると厳密に周期が求められない
ということです。なのでもし本当に求めたいなら、摩擦のない真空中で計測しましょう^^
振り子
円周率が求めたくなって、バネがない！そんな時でも そこに紐とボールさえがあれば、円周率を求めることができます！ 振り子のいいところは
- ばね定数などをあらかじめ測るべき定数がない．
というところ。バネはバネの種類によって周期が変わっちゃいますが、重力定数はほぼ普遍なので、どんなところでも使えます。
注意しないといけないのは、これは振り子の振れ幅が小さいという近似で成り立っているということ.
振り子の振れ幅を大きくしちゃうと、が成り立たなくなり、楕円関数を使わないといけないので注意しましょう！！
The Pi Machine
数年前、こんな論文が話題になりました
PLAYING POOL WITH π (THE NUMBER π FROM A BILLIARD POINT OF VIEW)
重さの違うボールをぶつけていくと、そのぶつかった回数が円周率になる。という論文です。
注意点は
- 完全弾性衝突のボールを用意する
- 精度良く質量比が求められている
- 空気抵抗がない環境を用意する
ことが必要です。これらの道具・環境が揃えられる人は是非やってみましょう！ 道具、環境を揃えるのが厳しい人は、シミュレーションでやってみましょう！
終わりに
いかがでしたか？単純に円周率、という以上に、様々な分野と深い関わりを見せていることがわかります。 たまにはこういうことに思いを馳せてみるのも楽しいですね！
*1:そういや、今日は国公立二次の入試試験の日ですね。受験生の方は、お疲れ様です。
