※ この記事は以下の動画の補足記事です。動画で触れられなかった内容、説明を飛ばした内容などを記載しています。
物体の衝突回数を数えると円周率が求められることが、2003年に発表された、東イリノイ大学の Gregory Galperin 氏による以下の論文で示されました。
PLAYING POOL WITH π (THE NUMBER π FROM A BILLIARD POINT OF VIEW)
この記事では、この論文に沿って円周率の計算方法を解説します。また、この論文では円周率の性質に関するとある予想 (おそらく未解決) が述べられているので、それについても解説します。
この方法では古典力学の基本的な事実である、運動量保存則と運動エネルギー保存則が用いられますが、記事の末尾に補足として、運動の三法則 (慣性の法則、作用・反作用の法則、運動方程式) $+ \alpha$ からそれらを導く方法を載せておきます。
計算のサンプルプログラムは以下で公開しています。
目次
概要
質量 $M$ kg の物体 $B_1$ と、質量 $m$ kg の物体 $B_2$ が一直線上に並んでいて、$B_1$ は $B_2$ と壁に挟まれているとします。
この状況で $B_1$ を壁の方向に動かすと、$B_1$ は $B_2$ に衝突し、$B_2$ は壁に向かって動きます。その後、$B_2$ は壁に当たって反発し、$B_1$ の方向に進みます。この後 $B_2$ が $B_1$ に衝突するかどうかは分かりませんが、条件によっては $B_2$ は壁および $B_1$ との衝突を何度か繰り返します。
$B_1$ と $B_2$ の質量比 $M / m$ が、ある自然数 $N$ が存在して $100^N$ と表されるとき、$B_2$ が壁および $B_1$ と衝突する回数は有限であり、衝突回数の数を数えると円周率の $N$ 桁目までの値が求まる、というのが論文の主定理です。
ただし、$B_1$, $B_2$ には (例えば摩擦や空気抵抗など、) 衝突による力以外の力は働かないとし、衝突は全て弾性衝突であるとします。壁との衝突時は、$B_1$ の速度の大きさはそのままで向きだけ変わるものとします。
以上の仮定のもとで、運動量は $B_1$ と $B_2$ の衝突で保存され、運動エネルギーは常に保存されるので (補足を参照)、$B_1$ の初速を $v$、ある時点での $B_1$ の速度を $v_1$, $B_2$ の速度を $v_2$ とおけば
\begin{gather} M v_1 + m v_2 = P \quad (\textrm{一定}) \\ \frac{1}{2} M v_1^2 + \frac{1}{2}m v_2^2 = \frac{1}{2} M V^2 = K \quad (\textrm{一定}) \\ \end{gather}
が成り立ちます (運動量は壁と $B_2$ との衝突では保存されない)。$X = \sqrt{m} v_1$, $Y = \sqrt{M} v_2$ とおけば
\begin{gather} \sqrt{m} X + \sqrt{M} Y = P \\ X^2 + Y^2 = 2K \end{gather}
と、直線と円の方程式になります。このように変形すると、円周率と関わりがあっても不思議では無いような気がします。
論文名に POOL (ビリヤード) と入っているようにこのモデルはビリヤードに例えられており、Galperin のビリヤードと呼ばれています。そのため元の論文では $B_1$, $B_2$ は球であるとされていますが、摩擦や空気抵抗などが無いので球である必要はありません。シミュレーションのアニメーション作る際、分かりやすいように重さの違いを物体の大きさの違いで表現しましたが、その場合は球にする衝突位置がずれてしまうのでブロック状にしました。
ちなみに量子力学版のモデルも考えられているようですが、この記事では解説しません。
Hear π from Quantum Galperin Billiards
衝突後の速度の計算
まずは、2物体の衝突後の速度が求められることを確認しましょう。$B_1$, $B_2$ の衝突前の速度をそれぞれ $v_1$, $v_2$、衝突前の速度を $v_1^{\prime}$, $v_2^{\prime}$ とおきます。このとき、運動量保存則と運動エネルギー保存則から
\begin{gather} M v_1 + m v_2 = M v_1^{\prime} + m v_2^{\prime} \\ \frac{1}{2}M {v_1}^2 + \frac{1}{2}m {v_2}^2 = \frac{1}{2}M {v_1^{\prime}}^2 + \frac{1}{2}m {v_2^{\prime}}^2 \end{gather}
が成り立ちます。運動エネルギー保存則の式を $m, M$ についてまとめると
\begin{align} M (v_1^2 -{v_1^{\prime}}^2) &= -m (v_2^2 -{v_2^{\prime}}^2) \\ M (v_1 -v_1^{\prime})(v_1 +v_1^{\prime}) &= -m (v_2 -v_2^{\prime})(v_2 +v_2^{\prime}) \end{align}
となりますが、運動量保存則の式を $M, m$ についてまとめると
$$M(v_1 -v_1^{\prime}) = m(v_2 -v_2^{\prime})$$
となり、衝突前後で速度が変わることから $v_1 \neq v_1^{\prime}$, $v_2 \neq v_2^{\prime}$ なので、両辺をそれぞれで割ると
$$v_1 +v_1^{\prime} = v_2 +v_2^{\prime}$$
が成り立つことがわかります。後はこの式と運動量保存の式から、$v_2^{\prime}$ を消せば $v_1^{\prime}$ が、$v_1^{\prime}$ を消せば $v_2^{\prime}$ が、$v_1$ と $v_2$ を用いて表されます。この記事では具体的な値はさほど重要ではありませんが、計算すると
\begin{align} v_1^{\prime} = \frac{(M -m) v_1 + 2 m v_2}{M+m} \\ v_2^{\prime} = \frac{ 2 M v_1 +(m -M) v_2}{M+m} \end{align}
となります。
配位空間
壁の位置が $0$、物体がある方向が正の方向となるように座標をとり、$B_1$ の座標を $x(t)$、$B_2$ の座標を $y(t)$ とします。簡単のため物体は大きさがないと仮定します。
そして組 $(x, y)$ をとり、2 物体の位置を $xy$ -平面 (配位空間という) 上の点とみなします。状況から $0 \leq y(t) \leq x(t)$ なので、2 物体の位置から定まる点は第1象限上の $y=x$ より下の部分にあります。点の進む方向は $(v_1, v_2)$ で与えられ、最初に $B_1$ を壁の方向に動かすと、点は左に動きます。
2 物体の衝突が起きるのは $y(t) = x(t)$ のとき、つまり点が斜めの直線にぶつかるときであり、その後点の進む方向が変わり、$(v_1, v_1)$ から $(v_1^{\prime}, v_2^{\prime})$ になります。$B_1$ が壁と衝突するのは $x = 0$ のとき、つまり点が $x$ 軸にぶつかるときであり、そのときに速度は $(v_1, v_2)$ から $(v_1, -v_2)$ になります。
このようにすると、直線上の 2 物体の運動が配位空間上の点の運動に置き換えられ、衝突回数は配位空間上の点と、直線 $y = x$ および $x$ 軸との接触回数と等しくなります。
ただし $y = x$ とぶつかった時の入射角と反射角は、一般には一致しません。例えば点の速度が $(-V, 0)$ の場合、つまり点が真っ直ぐ左に進む場合、入射角と反射角が一致するならば、点は衝突後真っ直ぐ下に進む、つまり $B_2$ の衝突後の速度 $v_2^{\prime}$ は $0$ となるはずですが、そのためには $M = m$ でなければいけません。この問題は、座標変換を行うことで解決します。
座標変換
$X = \sqrt{M}x$, $Y = \sqrt{m}y$ とおきます。つまり、$x$ 軸を $\sqrt{M}$ 倍拡大し、$y$ 軸を $\sqrt{m}$ 倍拡大します。すると、$y = x$ という曲線は
\begin{align} \frac{1}{\sqrt{m}} Y &= \frac{1}{\sqrt{M}} X \\ Y &= \frac{\sqrt{m}}{\sqrt{M}}X \end{align}
と変換されます。そして、$B_1$, $B_2$ の速度がそれぞれ $v_1$, $v_2$ であるとき、$XY$-座標における点の速度を $(V_1, V_2)$ とおくと
\begin{align} V_1 = \frac{d X}{dt} = \sqrt{M} \frac{dx}{dt} = \sqrt{M} v_1 \\ V_2 = \frac{d Y}{dt} = \sqrt{m} \frac{dy}{dt} = \sqrt{m} v_2 \end{align}
となります。これは軸を拡大しただけであることからもわかります。
ここで運動量保存則と運動エネルギー保存則を $XY$-座標で考えると、衝突前の速度を $V = (V_1, V_2)$, 衝突後の速度を $V^{\prime} = (V_1^{\prime}, V_2^{\prime})$ とおけば
\begin{gather} \sqrt{M} V_1 +\sqrt{m} V_2 = \sqrt{M} V_1^{\prime} +\sqrt{m} V_2^{\prime} \\ \frac{1}{2}{V_1}^{2} + \frac{1}{2}{V_2}^{2} = \frac{1}{2}{V_1^{\prime}}^{2} + \frac{1}{2}{V_2^{\prime}}^{2} \end{gather}
となります。これを内積を用いて表すと
\begin{gather} (\sqrt{M}, \sqrt{m}) \cdot V = (\sqrt{M}, \sqrt{m}) \cdot V^{\prime} \\ \frac{1}{2}||V||^2 = \frac{1}{2}||V^{\prime}||^2 \end{gather}
となります。
ベクトル $(\sqrt{M}, \sqrt{m})$ は $Y = \frac{\sqrt{m}}{\sqrt{M}} X$ に平行なベクトルであり、運動量保存則は $(\sqrt{M}, \sqrt{m})$ と速度の内積が不変である、運動エネルギー保存則は速度の大きさが不変であると言い換えられます。
ここで、上の図のように $V$ と $Y = \frac{\sqrt{m}}{\sqrt{M}} X$ の成す角を $\theta$、$V^{\prime}$ と $Y = \frac{\sqrt{m}}{\sqrt{M}} X$ の成す角を $\theta^{\prime}$ とおくと、運動量保存則から
\begin{align} & ||(\sqrt{M}, \sqrt{m})|| \cdot ||V|| \cos (\pi -\theta) \\ = \ & ||(\sqrt{M}, \sqrt{m})|| \cdot ||V^{\prime}|| \cos (\pi -\theta^{\prime})\end{align}
であり、運動エネルギー保存則から $||V|| = ||V^{\prime}||$ なので、$\cos (\pi -\theta) = \cos(\pi-\theta^{\prime})$ となります。$0 \leq \theta, \theta^{\prime} < \pi$ から $\theta = \theta^{\prime}$ がわかります。特に、入射角と反射角が等しいことがわかります。
$B_2$ と壁との衝突、$XY$-座標では点と $X$ 軸の衝突、では入射角と反射角が等しいことは明らかなので、どの衝突においても入射角と反射角は等しいです。
配位空間と衝突回数
後は、$XY$-平面において点の初期速度を $(-1, 0)$ としたとき、点が $Y = \frac{\sqrt{m}}{\sqrt{M}} X$ と $X$ 軸に何回ぶつかるのかを計算すれば良いです。ただし、素直に計算するのは難しいので少し工夫をします。
入射角と反射角が等しいことから、衝突後の軌跡を $Y = \frac{\sqrt{m}}{\sqrt{M}} X$ を対称軸として反転させたものは直線になります。
よって角度 $\alpha$ を
$$\tan \alpha = \frac{\sqrt{m}}{\sqrt{M}}$$
を満たすようにとったとき、衝突回数を $C_N$ とおけば
$$C_N \alpha < \pi \leq (C_N + 1) \alpha$$
を満たします。(下の図では $6$ 回衝突します。)
よって後は $\alpha$ の値を評価すれば、$\pi$ の値を $C_N$ を用いて近似することができます。
$\tan$ の逆関数は $\arctan$ と呼ばれており、
\begin{align} \tan (\arctan x) &= x \quad (x \in \mathbb{R}) \\ \arctan (\tan \theta) &= \theta \quad (-\frac{\pi}{2} < \theta < \frac{\pi}{2}) \end{align}
を満たします (正確には逆関数のうち、$(-\frac{\pi}{2}, \frac{\pi}{2})$ を値域とするもの) 。よって上の式の両辺の $\arctan$ をとれば
$$\alpha = \arctan \frac{\sqrt{m}}{\sqrt{M}} = \arctan (10^{-N})$$
となります。よって円周率は
$$C_N \arctan (10^{-N}) < \pi \leq (C_N + 1) \arctan (10^{-N}) \tag{1} $$
という不等式を満たし、$\arctan (10^{-N})$ の大きさを評価すれば円周率が求められます。
$\arctan$ の評価方法
積分表示
まずは $\arctan x$ の積分表示を与えます。$\theta = \arctan x$ とおくと $x =\tan \theta$ であり、逆関数の微分の公式から
\begin{align} (\arctan x)^{\prime} &= \frac{d \theta}{d x} = \frac{1}{\frac{d x}{d\theta}} \\ &= \frac{1}{(\tan \theta)^{\prime}} = \cos^2 \theta \\ &= \frac{\cos^2 \theta}{\sin^2 \theta + \cos^2 \theta} \\ &= \frac{1}{\tan^2 \theta + 1} = \frac{1}{x^2 + 1} \end{align}
となります。よって
$$\arctan x = \int_0^x \frac{1}{t^2 + 1} dt + C$$
となります。ここで $x = 0$ のとき、$0 = \tan \theta$ を満たす $\theta$ は $-\frac{\pi}{2} < \theta < \frac{\pi}{2}$ の範囲で $\theta = 0$ のみなので、$\arctan 0 = 0$ であり、$C = 0$ となります。従って
$$\arctan x = \int_0^x \frac{1}{t^2 + 1} dt$$
と $\arctan x$ を積分で表すことができます。
級数展開
$\frac{1}{t^2 + 1}$ を分解します。一般に
$$(1 -x)(1 +x +\cdots + {x^n}) = (1 +x^{n+1})$$
が成り立つので、$x \neq 1$ ならば両辺 $x -1$ で割って整理することで
$$\frac{1}{1 -x} = \sum_{k=0}^{n} x^k + \frac{x^{n+1}}{1 -x}$$
と表すことができます。$x = -t^2$ とおけば、$t^2 \neq -1$ なので
$$\frac{1}{1 +t^2} = \sum_{k=0}^{n} (-1)^k t^{2k} + (-1)^{n+1} \frac{t^{2n+2}}{1+ t^2}$$
となります。これをさっきの積分表示と合わせると
\begin{align} \arctan x &= \int_0^x \sum_{k=0}^{n} (-1)^k t^{2k} + (-1)^{n+1} \frac{t^{2n+2}}{1+ t^2} dt \\ &= \sum_{k=0}^{n} \frac{(-1)^k}{2k+1} x^{2k+1} + (-1)^{n+1} \int_0^x \frac{t^{2n+2}}{1+ t^2} dt \end{align}
となります。特に $n = 0$ とすれば
$$\arctan x = x -\int_0^x \frac{t^2}{1 +t^2} dt$$
となります。
ちなみにこれは
$$\frac{1}{1 +t^2} = 1 -\frac{t^2}{1+t^2}$$
を積分しても得られますし、さらに続けて
\begin{align} &\frac{1}{1 +t^2} = 1 -\frac{t^2}{1+t^2} \\ = \ & 1 -\frac{t^2 +t^4 -t^4}{1+t^2} = 1 -t^2 +\frac{t^4}{1+t^2}\\ = \ & 1 -t^2 +\frac{t^4 +t^6 -t^6}{1+t^2} = 1 -t^2 +t^4 -\frac{t^6}{1+t^2} \\ = \ & \cdots \end{align}
とすれば級数展開が得られます。
評価
$1 +t^2 > 1$ $(t > 0)$ であることから $x > 0$ のとき
$$0 < \int_0^x \frac{t^2}{1 +t^2} dt < \int_0^x t^2 dt = \frac{1}{3}x^3$$
が成り立ちます。よって
$$x -\frac{1}{3}x^3 < \arctan x < x \qquad (x > 0) \tag{2}$$
が成り立ちます。ここで $x = 10^{-N}$ とすれば
\begin{align} \frac{1}{10^N}\left(1 -\frac{1}{3 \cdot 10^{2N}}\right) < \arctan\frac{1}{10^N} < \frac{1}{10^N} \end{align}
が成り立ちます。衝突回数の不等式 $(1)$ と合わせれば
$$\frac{C_N}{10^N}\left(1 -\frac{1}{3 \cdot 10^{2N}}\right) < \pi < \frac{C_N}{10^N} +\frac{1}{10^N}$$
が成り立ちます。
円周率の $N$ 桁目までの値との関係
実数 $x$ に対して、$x$ 以下の整数のうち最大のものを $[x]$ と表すこととします。このとき衝突回数の不等式 $(1)$ から
$$C_N < \frac{\pi}{\arctan (10^{-N})} \leq C_N + 1$$
なので
$$C_N = \begin{cases}\left[ { \displaystyle \frac{\pi}{\arctan (10^{-N})} } \right] & \left({ \displaystyle \frac{\pi}{\arctan (10^{-N})}} \textrm{ は整数でない} \right) \\ {\displaystyle \frac{\pi}{\arctan (10^{-N})}} -1 & \left({\displaystyle \frac{\pi}{\arctan (10^{-N})}} \textrm{ は整数} \right) \end{cases}$$
が成り立ちます。
ここで、$\arctan x$ の不等式 $(2)$ から $0 < x \leq 1$ のとき
\begin{align} 0 &< \frac{1}{\arctan x} -\frac{1}{x} = \frac{x -\arctan x}{x \arctan x} \\ &< \frac{\frac{x^3}{3}}{x (x -\frac{x^3}{3})} = \frac{x}{3 -x^2} \leq \frac{x}{2} \end{align}
なので
$$\frac{1}{x} < \frac{1}{\arctan x} < \frac{1}{x} +\frac{x}{2} \qquad (0 < x \leq 1) \tag{3}$$
が成り立ちます。
まず ${ \displaystyle \frac{\pi}{\arctan (10^{-N})} }$ が整数のときは
\begin{align} & 0 < \frac{1}{\arctan x} -\frac{1}{x} \leq \frac{x}{2} \\ \Rightarrow \ & 0 < \frac{\pi}{\arctan (10^{-N})} -\frac{\pi}{10^{-N}} < \frac{\pi}{10^N} \end{align}
から、$N \geq 1$ ならば
$$[\pi \cdot 10^N] = { \displaystyle \frac{\pi}{\arctan (10^{-N})} } -1$$
なので
$$C_N = [\pi \cdot 10^N]$$
が成り立ちます。$N = 0$ ならば、$\arctan 1 = \frac{\pi}{4}$ であり、$C_N = 4 -1 = 3$ なので成り立ちます。
次に ${ \displaystyle \frac{\pi}{\arctan (10^{-N})} } $ が整数でないときは、不等式 $(3)$ から、$N \geq 1$ のとき
\begin{align} [\pi \cdot 10^N] \leq \left[ { \displaystyle \frac{\pi}{\arctan (10^{-N})} } \right] &\leq \left[\pi \cdot 10^N +\pi \cdot \frac{10^{-N}}{2} \right] \\ & \leq [\pi \cdot 10^N] + 1 \end{align}
となります ($N \geq 1$ のとき $\pi \cdot \frac{10^{-N}}{2} < 1$)。$N = 0$ のときは、$\arctan 1 = \frac{\pi}{4}$ なので $\frac{\pi}{\arctan 1}$ は整数になります。ここで、上の不等式からもし
$$\left[\pi \cdot 10^N +\pi \cdot \frac{10^{-N}}{2} \right] = [\pi \cdot 10^N]$$
が成り立てば
$$C_N = [\pi \cdot 10^N]$$
となります。上の条件が成り立たない、つまり
$$\left[\pi \cdot 10^N +\pi \cdot \frac{10^{-N}}{2} \right] = [\pi \cdot 10^N] + 1$$
が成り立つとすれば、$\pi$ を $10$ 進数展開したときの小数点以下 $N +1$ 桁目から $2N -1$ 桁目までが $9$ で、$2N$ 桁目が $8$ 以上である必要があります。特に小数点以下 $2N$ 桁目までに $9$ が $N-1$ 桁以上連続で続く必要があります。$N$ が大きくなるほど条件はキツくなるので、このようなことは起きないと思われます。つまり以下が成り立つと思われます。
予想 1.
$\pi$ を $10$ 進数で表したときの小数点以下 $n$ 桁目の値を $0 \leq a_n \leq 9$ とおく. このとき, 任意の整数 $N \geq 2$ に対して
$$a_{N+1} = a_{N+2} = \cdots = a_{2N-1} = 9$$
となることはない. $\Box$
予想 1 が正しいかどうかは未解決らしいですが、$N \leq 10^8$ までは正しいことが確認されているようです。そして、予想 1 が正しければ、以下の予想も正しいことがわかります ($N = 0, 1$ の場合は直接確かめれば良い)。
予想 2.
任意の整数 $N \geq 0$ に対して $C_N = [10^N \pi]$ が成り立つ。$\Box$
つまり、$C_N$ を求めることで $\pi$ の小数点以下 $N$ 桁目まで求めることができます。
補足（運動の三法則から保存則を導く）
運動の三法則
- 慣性の法則
- 作用・反作用の法則
- 運動方程式 $F = ma$
から運動量保存則と運動エネルギー保存則を導きます。
運動量保存則
$A_1$, $A_2$ を質点とし、この 2 質点の衝突を考えます。$A_1$ の質量、速度をそれぞれ $M$, $v_1$、$A_2$ の質量、速度をそれぞれ $m$, $v_2$ とします。このとき、各質点の運動量の和は
\begin{align} P &= M v_1 + mv_2 \end{align}
で与えられます。運動量の時間に関する微分は
\begin{align} \frac{dP}{dt} &= M \frac{d v_1}{dt} + m \frac{d v_2}{dt} \\ &= M a_1 + m a_2 \\ &= F_1 + F_2 \end{align}
となります ($F_i$ は $A_i$ に働く力で、時間 $t$ の関数) が、2 質点に働く力は 2 質点間の衝突によるものだけなので、作用・反作用の法則から
$$F_1 = -F_2$$
が成り立ちます。よって $\frac{dP}{dt} = 0$ となり、運動量 $P$ は保存されます。
運動エネルギー保存則
さっきと同様に、2 質点 $A_1$, $A_2$ の運動のみを考え、衝突前と衝突後で質点それぞれの運動エネルギーの和
$$K = \frac{1}{2}Mv_1^2 + \frac{1}{2}m v_2^2$$
が変わらないことを示します。ただし、$A_1$ と $A_2$ は完全弾性衝突するとします。
ここで完全弾性衝突するとは、衝突前後の相対速度の比が $-1$ であること、つまり
$$-1 = \frac{v_2^{\prime} -v_1^{\prime}}{v_2 -v_1}$$
を満たすことをいいます。相対的に見て、衝突前後で速度の大きさは変わらず向きだけ変わります。
衝突前の運動エネルギーの和を $K$、衝突後の運動エネルギーの和を $K^{\prime}$ とおくと
\begin{align} &K -K^{\prime} \\ = \ & \frac{1}{2}Mv_1^2 + \frac{1}{2}mv_2^2 -\left(\frac{1}{2}M {v_1^{\prime}}^2 + \frac{1}{2}m{v_2^{\prime}}^2\right) \\ = \ & \frac{1}{2}M(v_1^2 -{v_1^{\prime}}^2) +\frac{1}{2}m(v_2^2 -{v_2^{\prime}}^2) \\ = \ & \frac{1}{2}M(v_1-v_1^{\prime})(v_1+v_1^{\prime}) +\frac{1}{2}m(v_2 -v_2^{\prime})(v_2 +v_2^{\prime}) \end{align}
となりますが、完全弾性衝突の式を整理すると
$$v_1 + v_1^{\prime} = v_2 +v_2^{\prime}$$
となり、運動量保存則から
\begin{align} M v_1 + m v_2 &= M v_1^{\prime} +m v_2^{\prime} \\ \Rightarrow \ \ M (v_1 -v_1^{\prime}) &= -m(v_2 -v_2^{\prime}) \end{align}
なので、$K -K^{\prime} = 0$ つまり $K = K^{\prime}$ となります。よって衝突前後で運動エネルギーは保存されます。計算過程から、運動量保存則が成り立つという条件のもと、衝突前後で運動エネルギーが保存することと、その衝突が完全弾性衝突であることが同値であることがわかります。
参考文献
G. GALPERIN. PLAYING POOL WITH π (THE NUMBER π FROM A BILLIARD POINT OF VIEW)