Pythonプログラミング(ステップ4・反復処理・円周率の計算)

このページでは反復処理の例として円周率の計算について考える。

1.円周率の数値計算

円周率は、数学や自然科学・理工学のあらゆる領域に登場する、極めて重要な数であることは言うまでもなく、実際にそれを使った計算を行う機会が無くても、 π3.14 程度の値であることは誰でも知っている。

円周率は直径1の円の周長(半径1の円の周長の1/2)に他ならないから、円を多角形で近似した上で、正多角形の周の長さから見積もることができて、辺の数を増やすことで、近似の精度を上げることができる。 実際に、膨大な数の辺の正多角形を使った計算が続けられてきた。 こうした計算方法はとても手間がかかるため、一生をそれに費やした例さえもあったようである。

ニュートンの登場とともに、解析学的なアプローチによって、こうした計算は劇的に効率されることとなる。例えば、ニュートンは従来の計算方法を一新するような公式 π=334+24(232315251282717229570421171664213) を得ている。

ニュートンがゲームチェンジャーとなって、その後は、円周率を効率的に計算する公式や手法が開発されるようになり、 現在(2021年8月)では、スーパーコンピュータを使って、62兆8000億桁まで計算されているようである。 また、民生品のパソコンを使って5兆桁の計算を行なった例もある(こちらを参照)。

効率的な計算方法のうちで、比較的むかしから知られていなものに、マチンの公式 π=16arctan(15)4arctan(1239) がある。 ここで、逆正接関数の近似値はテイラー級数 arctan(x)n=0M(1)n2n+1x2n+1 によって四則演算だけで求めることができる。

decimalモジュール

Pythonには、任意の精度で計算を行うことができるdecimalモジュールが標準で備わっているので、それを用いてマチンの公式によって円周率を1000桁計算するコードの例の以下に示す。

# coding: utf-8
from decimal import *

getcontext().prec = 1000

def arctan(x,N):
    s = Decimal(0)
    p = x
    for n in range(N):
        k = Decimal(2*n+1)
        if n%2==0:
            s = s + Decimal(1)/k * p
        else:
            s = s - Decimal(1)/k * p
        p = p * x * x
    return s

m = 250
pi = Decimal(16)*arctan(Decimal(1)/Decimal(5), 3*m+2) - Decimal(4)*arctan(Decimal(1)/Decimal(239), m)

print(pi)

注: このコードで計算した値の最後の1桁は誤差を含んでいる。

icon-pc 練習:検算

数値計算で求めた円周率の値が正しいかどうかを検算するためにはどのような方法があるか、考えてみなさい。


2.非効率な計算方法

以下では、Pythonプログラミングの演習の題材として、敢えて効率の悪い円周率の計算方法を試してみることにしよう。


方法1.ライプニッツの公式

円周率と関係する無限級数は多数知られているが、比較的古くから知られている公式のひとつにライプニッツの級数 π4=n=0(1)n2n+1=113+1517+19 がある。 これは、等式 arctan(1)=π4 の左辺をテイラー級数で展開して得られる。

ニュートンの公式やマチンの公式とは違い、この交代級数はとても収束性が悪いことで知られている。

icon-pc 練習: ライプニッツ級数

ライプニッツの公式に基づいて円周率を計算するコードを書いて、動作を確認してみなさい。

icon-hint ヒント

上述のマチンの公式を使ったコードを再利用することができる。


方法2.写像の反復回数

漸化式 xn+1=xn2+c (ただし x0=0)を考える。

c>1/4xnは最終的に発散するが、十分に小さな ϵ>0を考えて、c=1/4+ϵと置くと、 値が「爆発」するまでに一定の猶予期間が生じることとなる。

そこで、あるϵに対して xn<2に留まっている状態で何回反復したかを N(ϵ)で表すことにすると、 N(ϵ)ϵπ となることが示せる。

すなわち、0から出発して、値が2を超えるまで反復した「回数」に ϵを乗じることで、円周率が得られる。

icon-pc 練習: 漸化式の発散までの回数

上記の方法で実際に反復回数を計算し、それに基づいてπを見積もるPythonコードを書いてみなさい。

icon-hint ヒント

参考文献:π in the Mandelbrot set, Aaron Klebanoff

関連動画(YouTube):Pi and the Mandelbrot Set - Numberphile


方法3.物体の衝突回数

高校で習った力学を思い出そう。ふたつの物体1と2が摩擦の無い床に置かれていて、x軸に沿って運動することができる。 x=0には壁があって、物体同士、そして物体と壁は完全弾性衝突(衝突の跳ね返り係数1)するとしよう。

ここで、物体1の質量をm1, 速度をv1と表記する。物体2についても同様に質量と速度を m2, v2 とする。

最初、物体2は静止しており(v2=0)、右側から物体1がそれに衝突する。そのときの速さをVとすれば、v1=Vである。

衝突される側(物体2)の質量が十分に小さければ(m2<m1)、物体2は左側に弾き飛ばされ、壁に衝突する。そして、壁で向きが180度反転して、 再び物体1と衝突する。

このように「小さな」物体2は、物体1と壁の間を何回か往復運動する過程で、物体1は右向きの運動量を得て、運動の向きが反転し、物体同士の衝突も終了する。

速度 v1,v2の物体同士が衝突する際、衝突後の速度 v1,v2は、運動量保存とエネルギー保存の式 m1v1+m2v2=m1v1+m2v212m1v12+12m2v22=12m1v12+12m2v22 を連立させることで解くことができる。

また、物体2が壁に衝突すると、単純に速度の向きが反転する、すなわち v2=v2 である。

n回目の衝突後の速度の組を(v1(n),v2(n))で表すことにすれば、以上の関係を使って、 (V,0)(v1(1),v2(1))(v1(2),v2(2)) のように速度が変化し、最終的に 0<v2<v1 となった状態で、もはや衝突は生じなくなる。

物体2の質量が十分に小さな場合(m2m1)、物体2の衝突回数(物体1との衝突と壁との衝突の総数)をNとすると、驚くべきことに Nm2m1π すなわち、衝突回数をカウントして、質量比の平方根を乗ずると、円周率の近似値が得られる。

icon-pc 練習: 衝突のシミュレーション

速度の変化を計算して衝突回数を見積り、円周率との関係性を確認してみなさい。

右の計算は、高校の物理のように反発係数 e を使って、 e=1=v1v2v1v2 とおいた上で、運動量保存の式と連立させたほうが簡単かもしれない。

icon-hint ヒント

参考動画(YouTube):The most unexpected answer to a counting puzzle, Why do colliding blocks compute pi?

衝突前後の速度変化の計算方法を以下におさらいしておこう。

系の運動量を P=m1v1+m2v2 エネルギーを E=12m1v12+12m2v22 と置く。 衝突後にもこれらの量が保存することから、衝突「後」の速度v1,v2は 連立方程式 m1v1+m2v2=P12m1v12+12m2v22=E を解いて得られる(この連立方程式が解を持たなければ、衝突は生じない)。

衝突が起こる場合は、簡単な計算から v2=bb24ac2a を得る。 ここで、 a=m2(m1+m2), b=2Pm2, c=P22Em1 と置いた。 これは、物体1との衝突によって物体2が壁に向かって動いている状態に相当する。 他方、物体1については v1=Pm2v2m1 である。

時間をおくと、物体2は(いずれ)壁に衝突するので、速度が反転して v2v2 となる。


方法4.モンテカルロ法

モンテカルロ法を使って、半径1の二次元球(円)の面積(すなわちπ)を求めることができる。 具体的な計算方法については、こちらのページを参照のこと。