「2つのボールをぶつけると円周率がわかる」のをしつこく確かめてみた・・・解析的に

前回の記事「「2つのボールをぶつけると円周率がわかる」らしいのでシミュレーションしてみた」の続き。　この記事では、衝突回数が円周率（の適当に桁をズラした値）になることを解析的に導いてみます。

• 大人になってからの再学習　『2つのボールをぶつけると円周率がわかる』

目次

• 問題設定 ▷
• n$n$ 回衝突後の速度 v(n)${\mathbf{\text{v}}}^{(n)}$ を求める（行列を用いた表式）▷
• Cm$C^m$ を求める▷
• v(n)${\mathbf{\text{v}}}^{(n)}$ を求める▷
• N$N$ と衝突回数 n$n$ との関係▷

問題設定

質量がそれぞれ 100Nm,m${100}^{N}m,m$ の2つの質点 m1,m2$m_1,m_2$ を考え、以下のような初期状態から運動が開始されるとします：



m1,m2$m_1,m_2$ の初速度はそれぞれ v,0$v,0$ とします（つまり質点2は最初静止している）。　質点同士の衝突と質点2と壁との衝突はどちらも完全弾性衝突（反発係数 e=1$e=1$）であるとし、衝突が起こらなくなるまでその衝突回数をカウントすると、その回数が [10Nπ]$[{10}^N\pi ]$ になることを示します。　ここで角括弧 [ ] はガウス記号*1です。

n$n$ 回衝突後の速度 v(n)${\mathbf{\text{v}}}^{(n)}$

回数をカウントする衝突には2種類あります：

• C1$C_1$ ： 質点同士の衝突
• C2$C_2$ ： 質点2と壁との衝突

前回見たように、この2種類の衝突は交互に起こります。　奇数回目の衝突は C1$C_1$、偶数回目の衝突は C2$C_2$ です。

n$n$ 回衝突した後の m1,m2$m_1,m_2$ の速度をそれぞれ v(n)1,v(n)2$v_1^{(n)},v_2^{(n)}$ とおきましょう（速度は2つの質点から壁の方向を正）：

これらを成分とした2成分のベクトルを考えます（要は速度空間を考える）：

　　v(n)=(v(n)1v(n)2)${\displaystyle \begin{array}{r}{\mathbf{\text{v}}}^{(n)}=\left(\begin{array}{c}{v}_1^{(n)}\\ v_2^{(n)}\end{array}\right)\end{array}}$

v(n)${\mathbf{\text{v}}}^{(n)}$ の漸化式と初期条件

前回の結果から、2つの衝突 C1,C2$C_1,C_2$ はこのベクトルに作用する行列として表すことができます。　奇数回目と偶数回目で衝突が異なるのでちょっと面倒ですが、C1,C2$C_1,C_2$ を行列と同一視して、以下のような v(n)$v^{(n)}$ の漸化式が得られます：

　　v(2m+1)v(2m)=C1v(2m)=C2v(2m−1)${\displaystyle \begin{array}{rl}{\mathbf{\text{v}}}^{(2m+1)}& =C_1{\mathbf{\text{v}}}^{(2m)}\\ {\mathbf{\text{v}}}^{(2m)}& =C_2{\mathbf{\text{v}}}^{(2m-1)}\end{array}}$

行列 C1,C2$C_1,C_2$ は前回計算していて

　　C1C2=(r−1r+12rr+12r+1−r−1r+1)=1r+1(r−12r2−r+1)=(100−1)${\displaystyle \begin{array}{rl}{C}_1& =\left(\begin{array}{cc}\frac{r-1}{r+1}& \frac{2}{r+1}\\ \frac{2r}{r+1}& -\frac{r-1}{r+1}\end{array}\right)=\frac{1}{r+1}\left(\begin{array}{cc}r-1& 2\\ 2r& -r+1\end{array}\right)\\ C_2& =\left(\begin{array}{cc}1& 0\\ 0& -1\end{array}\right)\end{array}}$

で与えられます。　ちなみに初期条件は

　　v0=(v0)${\displaystyle \begin{array}{r}{\mathbf{\text{v}}}_0=\left(\begin{array}{c}v\\ 0\end{array}\right)\end{array}}$

となります。

v(n)${\mathbf{\text{v}}}^{(n)}$ の一般項（行列を用いた表式）

行列 C$C$ を

　　C=C2C1=(r−1r+1−2rr+12r+1r−1r+1)=1r+1(r−1−2r2r−1)${\displaystyle \begin{array}{r}C=C_2{C}_1=\left(\begin{array}{cc}\frac{r-1}{r+1}& \frac{2}{r+1}\\ -\frac{2r}{r+1}& \frac{r-1}{r+1}\end{array}\right)=\frac{1}{r+1}\left(\begin{array}{cc}r-1& 2\\ -2r& r-1\end{array}\right)\end{array}}$

で定義すると、上記の初期条件と漸化式より

　　v(2m)v(2m+1)=C2C1C2C1C2C1⋯C2C1v0=(C2C1)mv0=Cmv0=C−12v(2m+2)=C2Cm+1v0(∵C−12=C2)${\displaystyle \begin{array}{rl}{\mathbf{\text{v}}}^{(2m)}& =C_2{C}_1{C}_2{C}_1{C}_2{C}_1\cdots C_2{C}_1{\mathbf{\text{v}}}_0\\ & =(C_2{C}_1{)}^m{\mathbf{\text{v}}}_0\\ & =C^m{\mathbf{\text{v}}}_0\\ {\mathbf{\text{v}}}^{(2m+1)}& =C_2^{-1}{\mathbf{\text{v}}}^{(2m+2)}\\ & =C_2{C}^{m+1}{\mathbf{\text{v}}}_0& (\because C_2^{-1}=C_2)\end{array}}$

となります。　v(2m+1)${\mathbf{\text{v}}}^{(2m+1)}$ は v(2m)${\mathbf{\text{v}}}^{(2m)}$ に C1$C_1$ を作用させても得られますが、上記の方が簡単。　よって、次に Cm$C^m$ を求めましょう。

Cm$C^m$ を求める

この節では Cm$C^m$ のもう少し簡単な形を求めます。

C1,C$C_1,C$ の書き換え

質点同士の衝突での速度の変換を表す行列 C1$C_1$ について、(1, 1) 成分と (2, 2) 成分が単に符号が異なるだけであることと、行列式が

　　|C1|=1(r+1)2∣∣∣r−12r2−r+1∣∣∣=1(r+1)2{−(r−1)2−4r}=−1${\displaystyle \begin{array}{rl}\left|C_1\right|& =\frac{1}{(r+1{)}^2}\left|\begin{array}{cc}r-1& 2\\ 2r& -r+1\end{array}\right|\\ & =\frac{1}{(r+1{)}^2}\{-(r-1{)}^2-4r\}\\ & =-1\end{array}}$

のように常に −1$-1$ になることを踏まえると、

　　sinθ=2r√r+1,cosθ=r−1r+1(∴tanθ=2r√r−1)${\displaystyle \begin{array}{rlrlr}\sin \theta & =\frac{2\sqrt{r}}{r+1},& \cos \theta & =\frac{r-1}{r+1}& (\therefore \tan \theta =\frac{2\sqrt{r}}{r-1})\end{array}}$

を満たす θ(0<θ≤π2)$\theta (0<\theta \le \frac{\pi }{2})$ と a=r√$a=\sqrt{r}$ を用いて、C1$C_1$ は

　　C1=(cosθasinθ1asinθ−cosθ)(tanθ=2r√r−1,a=r√)${\displaystyle \begin{array}{rlr}{C}_1& =\left(\begin{array}{cc}\cos \theta & \frac{1}{a}\sin \theta \\ a\sin \theta & -\cos \theta \end{array}\right)& (\tan \theta =\frac{2\sqrt{r}}{r-1},a=\sqrt{r})\end{array}}$

と書けます（r=1$r=1$ のときは θ=π2$\theta =\frac{\pi }{2}$）。　もともと C1$C_1$ は r$r$ だけで書けてたので、この θ$\theta$ と a$a$ は独立ではありませんが、それは後ほど。　C1$C_1$ をこの形に書いたとき、行列 C$C$ は

　　C=C2C1=(cosθ−asinθ1asinθcosθ)${\displaystyle \begin{array}{r}C=C_2{C}_1=\left(\begin{array}{cc}\cos \theta & \frac{1}{a}\sin \theta \\ -a\sin \theta & \cos \theta \end{array}\right)\end{array}}$

のように回転行列みたく書くことができます。

Cm$C^m$ を求める

上記の C$C$ の表式を用いると Cm$C^m$ が回転行列の冪乗のように簡単に計算できて

　　Cm=(cosmθ−asinmθ1asinmθcosmθ)${\displaystyle \begin{array}{r}{C}^m=\left(\begin{array}{cc}\cos m\theta & \frac{1}{a}\sin m\theta \\ -a\sin m\theta & \cos m\theta \end{array}\right)\end{array}}$

となります。　証明は数学的帰納法で。　m=1$m=1$ のときは成り立つのは明らか（m=0$m=0$ も単位行列になるので成り立ちますが）。　m=k$m=k$ のときに成り立つと仮定すると

　　Ck+1=CkC=(coskθ−asinkθ1asinkθcoskθ)(cosθ−asinθ1asinθcosθ)=(coskθcosθ−sinkθsinθ−asinkθcosθ−acoskθsinθ1acoskθsinθ+1asinkθcosθ−sinkθsinθ+coskθcosθ)=(cos(k+1)θ−asin(k+1)θ1asin(k+1)θcos(k+1)θ)${\displaystyle \begin{array}{rl}{C}^{k+1}& =C^{k}C\\ & =\left(\begin{array}{cc}\cos k\theta & \frac{1}{a}\sin k\theta \\ -a\sin k\theta & \cos k\theta \end{array}\right)\left(\begin{array}{cc}\cos \theta & \frac{1}{a}\sin \theta \\ -a\sin \theta & \cos \theta \end{array}\right)\\ & =\left(\begin{array}{cc}\cos k\theta \cos \theta -\sin k\theta \sin \theta & \frac{1}{a}\cos k\theta \sin \theta +\frac{1}{a}\sin k\theta \cos \theta \\ -a\sin k\theta \cos \theta -a\cos k\theta \sin \theta & -\sin k\theta \sin \theta +\cos k\theta \cos \theta \end{array}\right)\\ & =\left(\begin{array}{cc}\cos (k+1)\theta & \frac{1}{a}\sin (k+1)\theta \\ -a\sin (k+1)\theta & \cos (k+1)\theta \end{array}\right)\end{array}}$

よって m=k+1$m=k+1$ のときも成り立ちます。

v(n)${\mathbf{\text{v}}}^{(n)}$ を求める

以上をまとめると、n$n$ 回衝突した後の速度 v(n)${\mathbf{\text{v}}}^{(n)}$ が以下のように求まります：

　　v(2m)v(2m+1)=v(cosmθ−asinmθ)=v(cos(m+1)θasin(m+1)θ)${\displaystyle \begin{array}{rl}{\mathbf{\text{v}}}^{(2m)}& =v\left(\begin{array}{c}\cos m\theta \\ -a\sin m\theta \end{array}\right)\\ {\mathbf{\text{v}}}^{(2m+1)}& =v\left(\begin{array}{c}\cos (m+1)\theta \\ a\sin (m+1)\theta \end{array}\right)\end{array}}$

これで n$n$ 回衝突した後の速度が求まりました。　次は、N$N$ （これは r=100N$r={100}^N$ によって r$r$ と関係していて、a,θ$a,\theta$ は r$r$ を使って定義されていた）と衝突回数 n$n$ の関係を導きます。

N$N$ と衝突回数 n$n$ との関係

以下の手順によって、この関係を求めましょう：

1. 衝突が起こらなくなる条件を求める
2. ベクトル u(n)${\mathbf{\text{u}}}^{(n)}$ を導入
3. この条件から衝突回数 n$n$ と角度 α$\alpha$ （後述）の間の関係を求める
4. 衝突回数 n$n$ と N(r=100N)$N(r={100}^N)$ との関係を求める

衝突が起こらなくなる条件

この条件は前回見ました。

　　v1≤v2≤0${\displaystyle \begin{array}{r}{v}_1\le v_2\le 0\end{array}}$

衝突が n$n$ だけ起こるとき、v(n)${\mathbf{\text{v}}}^{(n)}$ が初めて上記の条件を満たせばいいことが分かります。

ベクトル u(n)${\mathbf{\text{u}}}^{(n)}$

ベクトル u(n)${\mathbf{\text{u}}}^{(n)}$ を

　　u(n)=(u(n)1u(n)2)=(v(n)11av(n)2)${\displaystyle \begin{array}{r}{\mathbf{\text{u}}}^{(n)}=\left(\begin{array}{c}{u}_1^{(n)}\\ u_2^{(n)}\end{array}\right)=\left(\begin{array}{c}{v}_1^{(n)}\\ \frac{1}{a}{v}_2^{(n)}\end{array}\right)\end{array}}$

で定義しましょう。　このとき

　　u(2m)u(2m+1)=v(cosmθ−sinmθ)=v(cos(m+1)θsin(m+1)θ)${\displaystyle \begin{array}{rl}{\mathbf{\text{u}}}^{(2m)}& =v\left(\begin{array}{c}\cos m\theta \\ -\sin m\theta \end{array}\right)\\ {\mathbf{\text{u}}}^{(2m+1)}& =v\left(\begin{array}{c}\cos (m+1)\theta \\ \sin (m+1)\theta \end{array}\right)\end{array}}$

となり、これらは u1u2$u_1{u}_2$-空間で円周

　　u21+u22=v2${\displaystyle \begin{array}{r}{u}_1^2+u_2^2=v^2\end{array}}$

上にあることが分かります。　また、衝突が起こらなくなる条件は

　　u1≤au2≤0${\displaystyle \begin{array}{r}{u}_1\le a{u}_2\le 0\end{array}}$

となります。　これらを図に描くと

となります。　角度 α$\alpha$ は a=r√$a=\sqrt{r}$ と

　　tanα=1a${\displaystyle \begin{array}{r}\tan \alpha =\frac{1}{a}\end{array}}$

で関係しています（0<α<π2$0<\alpha <\frac{\pi }{2}$）。　u(n)${\mathbf{\text{u}}}^{(n)}$ について、n$n$ が奇数のものを見ていくと、u(0)${\mathbf{\text{u}}}^{(0)}$ から円周の上半分を反時計回りに θ$\theta$ ずつ回転していきます。　また、n$n$ が偶数のものは、その1つ前の（n$n$ が奇数の） u(n)${\mathbf{\text{u}}}^{(n)}$ の真下に位置します。　つまり、円周の下半分を時計回りに θ$\theta$ ずつ回転していきます。　衝突が起こらなくなる条件は、この円周上の点が図で示した青色の領域（境界を含む）に入ることです。

補足

前節で θ,a$\theta ,a$ は r$r$ と

　　tanθ=2r√r−1a=r√${\displaystyle \begin{array}{rlr}\tan \theta & =\frac{2\sqrt{r}}{r-1}& a=\sqrt{r}\end{array}}$

と関係しているので θ,a$\theta ,a$ は互いに独立ではないと書きましたが、a$a$ の代わりに α$\alpha$ と θ$\theta$ との関係を見ておきましょう。　正接の倍角の公式（こちらを参照）より

　　tan2α=2tanα1−tan2α=2a1−1a2=2r√1−1r=2r√r−1=tanθ${\displaystyle \begin{array}{rl}\tan 2\alpha & =\frac{2\tan \alpha }{1-{\tan }^2\alpha }=\frac{2a}{1-\frac{1}{a^2}}=\frac{\frac{2}{\sqrt{r}}}{1-\frac{1}{r}}=\frac{2\sqrt{r}}{r-1}\\ & =\tan \theta \end{array}}$

となって、θ=2α$\theta =2\alpha$ であることがわかります。　したがって、行列 C$C$ は以下のように θ,α$\theta ,\alpha$ いずれかだけで書けます：

　　C=(cosθ−sinθtanθ2sinθcotθ2cosθ)=(cos2α−sin2αtanαsin2αcotαcos2α)${\displaystyle \begin{array}{r}C=\left(\begin{array}{cc}\cos \theta & \sin \theta \cot \frac{\theta }{2}\\ -\sin \theta \tan \frac{\theta }{2}& \cos \theta \end{array}\right)=\left(\begin{array}{cc}\cos 2\alpha & \sin 2\alpha \cot \alpha \\ -\sin 2\alpha \tan \alpha & \cos 2\alpha \end{array}\right)\end{array}}$

まぁ、この表式を使うわけではありませんが。

衝突回数 n$n$ と α$\alpha$ の関係を求める

まずは n=2m+1$n=2m+1$、つまり n$n$ が奇数の場合。　u(2m+1)${\mathbf{\text{u}}}^{(2m+1)}$ の偏角（u1$u_1$ 軸の正の部分から反時計回りに測った動径の角度）は (m+1)θ=2(m+1)α$(m+1)\theta =2(m+1)\alpha$ です。　このとき衝突が起こらなくなる条件は

　　∴π≤2(m+1)α<π+απα−1≤2m+1(=n)<πα${\displaystyle \begin{array}{rl}& \pi \le 2(m+1)\alpha <\pi +\alpha \\ \therefore & \frac{\pi }{\alpha }-1\le 2m+1(=n)<\frac{\pi }{\alpha }\end{array}}$

偏角が π+α$\pi +\alpha$ に等しくなる場合、2回前の 2m−1$2m-1$ 回衝突したときの偏角が π−α$\pi -\alpha$ となり、その次の 2m$2m$ 回衝突したときの偏角が π+α$\pi +\alpha$ となるので、この時点で衝突が終わります。

次は n=2m$n=2m$、つまり n$n$ が偶数の場合。　u(2m)${\mathbf{\text{u}}}^{(2m)}$ が円周の下半分にあることに注意して、これの偏角は 2π−mθ=2π−2mα$2\pi -m\theta =2\pi -2m\alpha$ です。　このとき衝突が起こらなくなる条件は

　　∴π<2π−2mα≤π+απα−1≤2m(=n)<πα${\displaystyle \begin{array}{rl}& \pi <2\pi -2m\alpha \le \pi +\alpha \\ \therefore & \frac{\pi }{\alpha }-1\le 2m(=n)<\frac{\pi }{\alpha }\end{array}}$

偏角が π$\pi$ に等しくなる場合、その1回前の 2m$2m$ の衝突の時点で偏角が π$\pi$ になるので、それ以上衝突は起こりません。

以上をまとめると、衝突回数が偶数、奇数どちらの場合もまとめて書けて

　　πα−1≤n<πα${\displaystyle \begin{array}{r}\frac{\pi }{\alpha }-1\le n<\frac{\pi }{\alpha }\end{array}}$

もしくはガウス記号 [ ] を用いて

　　n=[πα]${\displaystyle \begin{array}{r}n=\left[\frac{\pi }{\alpha }\right]\end{array}}$

となります。

衝突回数 n$n$ と N$N$ との関係を求める

α$\alpha$ が充分小さいとき（N$N$ が充分大きいとき）、

　　α∼tanα=1a=1r√=110N${\displaystyle \begin{array}{r}\alpha \sim \tan \alpha =\frac{1}{a}=\frac{1}{\sqrt{r}}=\frac{1}{{10}^N}\end{array}}$

が成り立つので、結局衝突回数 n$n$ は

　　n=[10Nπ]${\displaystyle \begin{array}{r}n=\left[{10}^N\pi \right]\end{array}}$

で与えられます。　前回行ったシミュレーションでは、N$N$ が充分大きいという条件は必要なく、0以上の自然数で成り立っていましたが、評価が面倒なのでこの辺で終了。

まとめ

思ったより長い導出になりましたが（最後の近似の評価は投げたけど）、一応高校の数学と物理の知識で大体出せましたね。　と言いつつ、行列は高校の数学からなくなったらしいですが。　まぁこういう話は、導出はある程度理解できても、そもそも「衝突回数が円周率を与える」ということに気付くなんてことがそう簡単にできねぇよ、という結論になってしまうんですが。　原論文とか読めばそのあたりのことも書かれてるのかもしれませんが、あんまり読む気しない・・・　ビリヤードっぽい基本的な力学系なら、量子ビットか何かの関連で量子系への拡張することを念頭においてまずは古典系でやってみた、みたいなもんじゃないかと推測してみたり。

ちなみに

この記事は海外では1ヶ月前くらいに雑誌に掲載されたとか書いてましたが、3月といえば3月14日の円周率の日にネタとして載ってたのかな？

不思議な数πの伝記

*1:挟まれた数を超えない最大の整数を与える。