プロクラシスト

今日の寄り道 明日の近道

もう円周率で悩まない!πの求め方10選


スポンサーリンク

こんにちは!ほけきよです。

皆さん、πを知っていますか??あの3.14以降無限に続く円周率です。

昔、どこかのお偉いさんが「3.14って中途半端じゃね?www3にしようぜ」

とかいって一時期円周率が3になりかけました。でもそれは円じゃなくて六角形だからだめです。全然ダメ。

それを受けて「あほか、円周率をちゃんと教えろ」 と主張したのが東大のこの問題*1

f:id:imslotter:20161119092147p:plain

めっちゃ単純な問題。でも、東大受験生でさえ 「普段強制的に覚えさせられたπというやつ、どうやったら求められるの???」 と悩んだことでしょう。

また、普段生活してると

「π求めてぇ」

と悩むこともあるでしょう。今日はそんなみなさんに、様々なπの求め方をお教えします。これで、あらゆる状況で求められるようになりますよ!

東大の問題へのアプローチ2つ

もちろん、πの厳密な値を求めることはできません。今でもπの値は日々計算され続けています。 じゃあ、πより少し小さい値で、うまくπの値を近似できる方法を考えよう。 というアプローチです。

多角形で近似

おそらく一番多かったであろう回答が、この多角形近似です

同じ半径であれば、正多角形はすべて円の中に収まります。正方形も正六角形も正八角形も。 なので、それを利用してやりましょう。正六角形は周と直径の比が3であることは簡単にわかるので

  • 正六角形よりも多角形
  • sinやcosの値が出せそう

な正八角形(もしくは正十二角形)を選びます。

解法はこんな感じです。

f:id:imslotter:20170225151832p:plain

tanの逆関数を使う

この問題に関しては、こんな解法もできます! 高3のときに習いますね!

  • 置換積分を使うと、答えにπが現れる
  • かつ、上に凸な関数
  • かつ、値を代入した時に計算がしやすい

と言えば、そう、 {
f(x) = \frac{1}{1+x^{2}}
}ですね!!

{
f(x) = \frac{1}{\sqrt{1-x^{2}}}
}は、ルートがある分、ちと使いにくいのです。 解法は↓のような感じ f:id:imslotter:20170225152746p:plain

無限級数を覚えておく

フーリエ級数を用いる

世の中にはこんな不思議な式があります

{
\sum _{{n=1}}^{{\infty }}{\frac  {1}{n^{{2}}}}=\frac{1}{1^{2}}+\frac{1}{2^{2}}+\frac{1}{3^{2}}+\cdots={\frac  {\pi ^{2}}{6}}
}

これを理解するためには,Fourier級数を知る必要があります。理系の方なら大学1-2年くらいで学びますね。 打ち切り項数と{\pi}の関係はこんな感じ。

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

フーリエ級数がわかれば、上の式以外にも、例えばこんな式も作れるようになります

{
\sum _{n=0}^{\infty }{\frac {(-1)^{n}}{2n+1}}=1-\frac{1}{3}+\frac{1}{5}-\frac{1}{7}+\cdots={\frac {\pi }{4}}
}

分数なら簡単に計算できるし,πも簡単に求められそうですね^^

ラマヌジャン式を使う

無性にπが求めたくなった時も,この無限級数を知っているだけでOK! あの天才ラマヌジャンが導出した式です

{
\frac{4}{\pi}=\sum _{{n=0}}^{\infty}{\frac{(-1)^{n}(4n)!(1123+21460n)}{882^{2n+1}(4^{n}n!)^{4}}}
}

美しい式ですね(白目)

めちゃくちゃ収束が早いことが知られているので,n=0,1,2とかをぶち込んでやるだけでそれなりの精度が出るのがいいところ

n = 0, 1での代入結果がこちら

n:0     Value:3.14158504007123751123
n:1     Value:3.14159265359762196468

n=0で、もう良さげ。すごい精度。

  • ちょっと複雑で覚えにくい
  • 分子分母の値がでっかくなりすぎて計算がそもそも厳しい

のがたまに傷かな??

コンピュータを使う

モンテカルロサンプリングする

あなたの眼の前にそこそこいいパソコンがあるなら,モンテカルロサンプリングでπを求めましょう!

  • {(x,y)} を,両方とも-1から1の範囲でランダムに選択
  • {x^{2}+y^{2}}を計算
  • {x^{2} + y^{2} \leq 1} なら +1, {x^{2}+y^{2} \gt 1} なら何もしない
  • N回繰り返して点をばらまく
  • {x^{2}+y^{2} \leq 1}だった点の数をNで割る

最終的にこの結果を4倍すればPiが求められます

いいところは,回数をこなせばこなすほど精度が上がるところと、事前に初期値設定が必要ないところ。

f:id:imslotter:20170226120649p:plain
点を打つほど円がわかりやすくなってくる

悪いところは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桁(約2兆6000億桁)の計算がされた(Wikipediaより)

なんかすごそう…よっぽど複雑なのかと思いきや、アルゴリズムは超簡単(Wikipediaより)

f:id:imslotter:20170225164803p:plain

実際にコードを書いてみて動かした結果がこちら

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回の更新でモンテカルロサンプリングを超えていることがわかります。しかも更新も一瞬! かなり優秀なアルゴリズムのようです。

実験で求める

ビュフォンの針

もしあなたが針やつまようじを大量に持っているならば、こんな実験をしてみましょう

  • 針の長さを測る({h})
  • その半分の長さの間隔で平行線を引く({\frac{h}{2}})
  • 引いた平行線の上に針をぶちまける
  • 平行線に交わる針の数を数える

f:id:imslotter:20170226130452p:plain

これはビュフォンの針問題と言って、針の数をめちゃくちゃ増やすと {p=\frac{1}{\pi}}となります。

こうするだけで、なんと{\pi}が求まります。ね、簡単でしょ???

単振動

円周率が求めたいときに、バネを見つけたとします。 それはラッキーですね。早速バネの振動する周期を求めましょう!!

図のように、周期に{\pi}が含まれているので、ばねの振動する時間を求めるだけで、簡単に{\pi}が求まります。 f:id:imslotter:20170226131947p:plain

注意点は

  • 摩擦があると厳密に周期が求められない
  • 空気抵抗があると厳密に周期が求められない

ということです。なのでもし本当に求めたいなら、摩擦のない真空中で計測しましょう^^

振り子

円周率が求めたくなって、バネがない!そんな時でも そこに紐とボールさえがあれば、円周率を求めることができます! 振り子のいいところは

  • ばね定数などをあらかじめ測るべき定数がない.

というところ。バネはバネの種類によって周期が変わっちゃいますが、重力定数はほぼ普遍なので、どんなところでも使えます。

f:id:imslotter:20170226132000p:plain

注意しないといけないのは、これは振り子の振れ幅が小さいという近似で成り立っているということ. 振り子の振れ幅を大きくしちゃうと、{\sin{\theta} \approx \theta }が成り立たなくなり、楕円関数を使わないといけないので注意しましょう!!

The Pi Machine

数年前、こんな論文が話題になりました

PLAYING POOL WITH π (THE NUMBER π FROM A BILLIARD POINT OF VIEW)

重さの違うボールをぶつけていくと、そのぶつかった回数が円周率になる。という論文です。

www.youtube.com

注意点は

  • 完全弾性衝突のボールを用意する
  • 精度良く質量比が求められている
  • 空気抵抗がない環境を用意する

ことが必要です。これらの道具・環境が揃えられる人は是非やってみましょう! 道具、環境を揃えるのが厳しい人は、シミュレーションでやってみましょう!

終わりに

いかがでしたか?単純に円周率、という以上に、様々な分野と深い関わりを見せていることがわかります。 たまにはこういうことに思いを馳せてみるのも楽しいですね!

魅惑のπ。
他に面白い求め方を知っている人は、教えてください!ではでは!

*1:そういや、今日は国公立二次の入試試験の日ですね。受験生の方は、お疲れ様です。

PROCRASIST