Pythonプログラミング(ステップ4・反復処理・円周率の計算)
このページでは反復処理の例として円周率の計算について考える。
1.円周率の数値計算
円周率は、数学や自然科学・理工学のあらゆる領域に登場する、極めて重要な数であることは言うまでもなく、実際にそれを使った計算を行う機会が無くても、
円周率は直径1の円の周長(半径1の円の周長の1/2)に他ならないから、円を多角形で近似した上で、正多角形の周の長さから見積もることができて、辺の数を増やすことで、近似の精度を上げることができる。 実際に、膨大な数の辺の正多角形を使った計算が続けられてきた。 こうした計算方法はとても手間がかかるため、一生をそれに費やした例さえもあったようである。
ニュートンの登場とともに、解析学的なアプローチによって、こうした計算は劇的に効率されることとなる。例えば、ニュートンは従来の計算方法を一新するような公式
ニュートンがゲームチェンジャーとなって、その後は、円周率を効率的に計算する公式や手法が開発されるようになり、 現在(2021年8月)では、スーパーコンピュータを使って、62兆8000億桁まで計算されているようである。 また、民生品のパソコンを使って5兆桁の計算を行なった例もある(こちらを参照)。
効率的な計算方法のうちで、比較的むかしから知られていなものに、マチンの公式
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桁は誤差を含んでいる。
練習:検算
数値計算で求めた円周率の値が正しいかどうかを検算するためにはどのような方法があるか、考えてみなさい。
2.非効率な計算方法
以下では、Pythonプログラミングの演習の題材として、敢えて効率の悪い円周率の計算方法を試してみることにしよう。
方法1.ライプニッツの公式
円周率と関係する無限級数は多数知られているが、比較的古くから知られている公式のひとつにライプニッツの級数
ニュートンの公式やマチンの公式とは違い、この交代級数はとても収束性が悪いことで知られている。
練習: ライプニッツ級数
ライプニッツの公式に基づいて円周率を計算するコードを書いて、動作を確認してみなさい。
ヒント
上述のマチンの公式を使ったコードを再利用することができる。
方法2.写像の反復回数
漸化式
そこで、ある
すなわち、0から出発して、値が2を超えるまで反復した「回数」に
練習: 漸化式の発散までの回数
上記の方法で実際に反復回数を計算し、それに基づいて
ヒント
参考文献:π in the Mandelbrot set, Aaron Klebanoff
関連動画(YouTube):Pi and the Mandelbrot Set - Numberphile
方法3.物体の衝突回数
高校で習った力学を思い出そう。ふたつの物体1と2が摩擦の無い床に置かれていて、x軸に沿って運動することができる。
ここで、物体1の質量を
最初、物体2は静止しており(
衝突される側(物体2)の質量が十分に小さければ(
このように「小さな」物体2は、物体1と壁の間を何回か往復運動する過程で、物体1は右向きの運動量を得て、運動の向きが反転し、物体同士の衝突も終了する。
速度
また、物体2が壁に衝突すると、単純に速度の向きが反転する、すなわち
物体2の質量が十分に小さな場合(
練習: 衝突のシミュレーション
速度の変化を計算して衝突回数を見積り、円周率との関係性を確認してみなさい。
右の計算は、高校の物理のように反発係数
ヒント
参考動画(YouTube):The most unexpected answer to a counting puzzle, Why do colliding blocks compute pi?
衝突前後の速度変化の計算方法を以下におさらいしておこう。
系の運動量を
衝突が起こる場合は、簡単な計算から
時間をおくと、物体2は(いずれ)壁に衝突するので、速度が反転して
方法4.モンテカルロ法
モンテカルロ法を使って、半径1の二次元球(円)の面積(すなわち