「2つのボールをぶつけると円周率がわかる」らしいのでシミュレーションしてみた

「2つのボールをぶつけると円周率がわかる」

だそうなので、ちょっと便乗記事を。

衝突回数が円周率（の桁をズラした値）になることの数学的証明まで追うのは大変そうなので、この記事では

衝突による速度の変化と衝突しなくなる条件を数式で書き下す
衝突回数はシミュレーションで求める

という方針で。

衝突による速度の変化

回数を数える衝突には以下の2つがあります：

質点同士の衝突
質点と壁の衝突

それぞれがどのように速度を変化させるのかを見ていきましょう。　2つの質点から壁に向かう方向を速度の正の方向とします。

質点同士の衝突

参照記事には書いてませんが、2つの質点の衝突は完全弾性衝突（反発係数 ${ e = 1 }$ ）としていいんだと思います。　最初に動かす方の質点を ${ m_1 }$ 、最初は止まっている方の質点を ${ m_2 }$ とし、それぞれの質量を ${ rm,\, m \; (r = 100^N) }$ とします（このとき衝突回数が ${ [10^N\pi] }$ *1 になるというのが元記事で証明されたことだそう）。　これらの質点がそれぞれ速度 ${ v_1,\,v_2 }$ で運動していて、完全弾性衝突によって速度が変化した後に ${ v_1',\, v_2' }$ となったとすると、運動量保存則より

　　 ${ \displaystyle\begin{align*} rmv_1' + mv_2' = rmv_1 + mv_2 \end{align*}}$

また、完全弾性衝突という条件より

　　 ${ \displaystyle\begin{align*} -\frac{v_2' - v_1'}{v_2 - v_1} = 1 \end{align*}}$

これらを ${ v_1', \, v_2' }$ を変数として連立方程式を解くと

　　 ${ \displaystyle\begin{align*} \begin{cases} v_1' &= \frac{r-1}{r+1}v_1 + \frac{2}{r+1}v_2 \\ v_2' &= \frac{2r}{r+1}v_1 - \frac{r-1}{r+1}v_2 \end{cases} \end{align*}}$

この記事では使いませんが、速度の変換を行列で書くと

　　 ${ \displaystyle\begin{align*} \begin{pmatrix} v'_1 \\ v'_2 \end{pmatrix} = \frac{1}{r+1} \begin{pmatrix} r-1 & 2 \\ 2r & -r+1 \end{pmatrix} \begin{pmatrix} v_1 \\ v_2 \end{pmatrix} \end{align*}}$

となります。

ちなみに、この衝突が起こるためには（2つの質点の位置関係は変わらないので）質点1の速度は質点2の速度より大きくなくてはなりません。　つまり ${ v_1 - v_2 > 0 }$ 。　このとき、反発係数の式より

　　 ${ \displaystyle\begin{align*} v_2' - v_1' = - (v_2 - v_1) > 0 \end{align*}}$

となって、今度は質点2の速度は質点1の速度より大きくなります。　したがって、2つの質点同士が衝突した後は、（もし衝突が起こるなら）次の衝突は質点2と壁との衝突になります。

質点と壁の衝突

壁と衝突する質点は質点2だけです。　壁との衝突も完全弾性衝突（反発係数 ${ e = 1 }$ ）として速度 ${ v_2 }$ で運動していたものが速度 ${ v_2' }$ になったとすると

　　 ${ \displaystyle\begin{align*} v_2' = -v_2 \end{align*}}$

これも行列の形で書いておきましょう（質点1の速度が変わらないのも考慮して）

　　 ${ \displaystyle\begin{align*} \begin{pmatrix} 1 & 0 \\ 0 & -1 \end{pmatrix} \end{align*}}$

まぁ、やっぱり使いませんが。

衝突しなくなる条件

2つの質点同士や質点2と壁との衝突が起こらなくなる条件は

質点1、質点2の速度がどちらも負
質点1の速度が質点2の速度よりも小さい（速度の絶対値、すなわち速さで比べると、質点1の方が大きい）

式で書くと

　　 ${ \displaystyle\begin{align*} v_1 \le v_2 \le 0 \end{align*}}$

質点2が止まった場合や、質点1と質点2の速度が（負の値で）等しいときも、もう衝突が起こりませんね。

Groovy コード

さて、以上の条件から衝突回数を解析に解いて出せればいいんですが、ここでは衝突をシミュレーションしてみましょう。　アニメーションは無しで。　質点同士の衝突後の速度を計算する部分が少々複雑に見える以外は結構簡単に書けます。　コードは Groovy で（Java, Scala, Go, Rust についてはこちら）。

def n = args.size() > 0 && args[0].isInteger() ? args[0].toInteger() : 2
def r = 100**n    // 質点1の質点2に対する質量比 m1 / m2 = 100^n

// 質点同士の衝突
def collide1 = { v -> [(r-1)*v[0] + 2*v[1], 2*r*v[0] - (r-1)*v[1]].collect{ it/(r+1) } }

// 質点2と壁との衝突
def collide2 = { v -> [v[0], -v[1]] }

def v = [1, 0]    // (v1, v2) で、初速度は (v1, v2) = (1, 0)
def count = 0
def collide = [collide1, collide2]

println "N = $n"
//println "count : [v1, v2]"
//println "$count : $v"

while(!stop(v)){
    v = collide[0](v)
    count++
    collide = collide.reverse()    // 質点同士の衝突と壁との衝突を交互に
//    println "$count : $v"
}

println "The number of collisions is $count."

def boolean stop(v){
    return v[0] <= v[1] && v[1] <= 0
}

引数に ${ N }$ の値（質点1と質点2の比が ${ m_1 : m_2 = 100^N : 1 }$ になる）を指定して実行します（デフォルト値は2）。　シミュレート中の速度の値を書き出したい場合は「println」文の前についているコメントを外して実行して下さい。　いくつかの ${ N }$ の値で試すと（上記のコードのソースは pi.groovy とします）

>groovy pi 0
N = 0
The number of collisions is 3.

>groovy pi 1
N = 1
The number of collisions is 31.

>groovy pi 2
N = 2
The number of collisions is 314.

>groovy pi 3
N = 3
The number of collisions is 3141.

>groovy pi 4
N = 4
The number of collisions is 31415.

>groovy pi 5
N = 5
The number of collisions is 314159.

おぉ、ちゃんと円周率の最初の数字になっとる。　他の言語でも練習に書いてみたいレベルのコードやね。

ちなみに「「2つのボールをぶつけると円周率がわかる」のをしつこく確かめてみた・・・解析的に」にてもうちょっと手で解いてみました。