こんにちは！ほけきよです。
皆さん、πを知っていますか？？あの3.14以降無限に続く円周率 です。
昔、どこかのお偉いさんが「3.14って中途半端じゃね？ｗｗｗ3にしようぜ」
とかいって一時期円周率が3になりかけました。でもそれは円じゃなくて六角形 だからだめです。全然ダメ。
それを受けて「あほか、円周率をちゃんと教えろ」
と主張したのが東大のこの問題*1
めっちゃ単純な問題。でも、東大受験生でさえ
「普段強制的に覚えさせられたπというやつ、どうやったら求められるの???」
と悩んだことでしょう。
また、普段生活してると
「π求めてぇ」
と悩むこともあるでしょう。今日はそんなみなさんに、様々なπの求め方をお教えします。これで、あらゆる状況で求められるようになり ますよ！
東大の問題へのアプローチ2つ
もちろん、πの厳密な値を求めることはできません。今でもπの値は日々計算され続けています。
じゃあ、πより少し小さい値で、うまくπの値を近似できる方法を考えよう。
というアプローチです。
多角形で近似
おそらく一番多かったであろう回答が、この多角形近似
同じ半径であれば、正多角形はすべて円の中に収まります。正方形も正六角形も正八角 形も。
なので、それを利用してやりましょう。正六角形は周と直径の比が3であることは簡単にわかるので
正六角形よりも多角形
sinやcosの値が出せそう
な正八角 形(もしくは正十二角形)を選びます。
解法はこんな感じです。
tanの逆関数 を使う
この問題に関しては、こんな解法もできます！
高３のときに習いますね！
置換積分 を使うと、答えにπが現れる
かつ、上に凸な関数
かつ、値を代入した時に計算がしやすい
と言えば、そう、
無限級数 を覚えておく
世の中にはこんな不思議な式があります
これを理解するためには，Fourier級数
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で、もう良さげ。すごい精度。
ちょっと複雑で覚えにくい
分子分母の値がでっかくなりすぎて計算がそもそも厳しい
のがたまに傷かな？？
コンピュータを使う
あなたの眼の前にそこそこいいパソコンがあるなら，モンテカルロ サンプリング
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
無限に時間と計算資源がある人は，試してみましょう！
もっと精度よく効率的に求めたい！！というアナタ！
ガウス =ルジャンドル のアルゴリズム
ガウス＝ルジャンドルのアルゴリズム - Wikipedia
ガウス ＝ルジャンドル のアルゴリズム は円周率を計算する際に用いられる数学の反復計算アルゴリズム である。円周率を計算するものの中では非常に収束が速く、2009年にこの式を用いて2,576,980,370,000桁 Wikipedia より)
なんかすごそう…よっぽど複雑なのかと思いきや、アルゴリズム は超簡単(Wikipedia より)
実際にコードを書いてみて動かした結果がこちら
import numpy as np
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
a = 1.0
b = 1 /np.sqrt(2 )
t = 0.25
p = 1.0
print ("0 : {0:.10f}" .format ((a+b)**2 /(4 *t)))
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)
重さの違うボール をぶつけていくと、そのぶつかった回数が円周率になる
VIDEO www.youtube.com
注意点は
完全弾性衝突のボールを用意する
精度良く質量比が求められている
空気抵抗がない環境を用意する
ことが必要です。これらの道具・環境が揃えられる人は是非やってみましょう！
道具、環境を揃えるのが厳しい人は、シミュレーション でやってみましょう！
終わりに
いかがでしたか？単純に円周率、という以上に、様々な分野と深い関わりを見せていることがわかります。
たまにはこういうことに思いを馳せてみるのも楽しいですね！
魅惑のπ。
他に面白い求め方を知っている人は、教えてください！ではでは！
