弾性衝突で円周率を求める話

最近、YouTubeでこんな動画を見かけました。

2つの物体を衝突させる話です。左端には壁があり、2つの物体を並べ、右側の物体を左に向かって滑らせます。摩擦はなく、すべての衝突が弾性衝突だとすると、左と右の物体の質量比が$100^n$のとき、左の物体の衝突回数は円周率の小数第$n$位までを整数で表したものになるというのです。

少し調べてみると比較的有名なことのようですが、私は初耳でした。私はうっかり解説動画から見てしまったので、自分で考える間もなくなぜそうなるかを知ってしまったのですが、もしも初めて聞いた人は答えを見ずに自分で考えてみることをお勧めします。何もヒントなしに証明できた人がいたとしたら大したものです。

自分で考えてみる 

さて、解説動画を見ているとところどころで天才的なひらめきが出てくるのですが、いったんそれは忘れて少し考えてみましょう。

右向きに$x$軸を取り、右側の物体を物体1とし、その質量と速度をそれぞれ$m_1$, $v_1$、左側の物体を物体2としｋ、その質量と速度をそれぞれ$m_2$, $v_2$と置きます。

こういうのはそれぞれの物体の位置と速度を計算するとドツボにはまります。使用すべきはエネルギー保存則と運動量保存則です。2つの物体が衝突するときは弾性衝突なので2つの物体のエネルギーの合計も運動量の合計も保存します。左の物体が壁にぶつかるときは、弾性衝突なのでエネルギーは保存しますが運動量は正負反転します。

まずは2つの物体が衝突するときのことを考えてみましょう。

\[\left\{
\begin{array}{l}
\dfrac{1}{2}m_1v_1^2+\dfrac{1}{2}m_2v_2^2=k \\
m_1v_1+m_2v_2=p
\end{array}
\right. \]

何の変哲もない、ただのエネルギー保存則と運動量保存則です。私は凡人なので、これを連立させて2次方程式を解こうとしました。$v_2$を削除し$v_1$を求めます。ただの2次方程式なので解の公式を使えば求められますが、結構計算量が多くしんどい計算となります。

\[v_1=\dfrac{p\pm{m_2\sqrt{2k(\frac{1}{m_1}+\frac{1}{m_2})-\frac{1}{m_1m_2}p^2}}}{m_1+m_2}\]

ルートの中身を$m_1$と$m_2$が対称になるように整理するのがポイントです。これを運動量保存則の式に当てはめて$v_2$を計算します。

\[v_2=\dfrac{p\mp{m_1\sqrt{2k(\frac{1}{m_1}+\frac{1}{m_2})-\frac{1}{m_1m_2}p^2}}}{m_1+m_2}\]

ルートの中身を対称にしたおかげでかなりきれいに計算できました。

この式が意味するところは、「2つの物体のエネルギーの合計が$k$、運動量の合計が$p$の場合、$v_1$と$v_2$の取りうる組み合わせは2ペアある」ということになります。これはすなわち物体の衝突前の状態と衝突後の状態です。衝突の後は物体1はより右向きの速度が上がりますので、下側の符号が衝突前、上側の符号が衝突後の速度ということになりますね。

さて、物体同士が衝突した後に物体2が壁にぶつかると物体2の運動量の正負が反転します。すなわち、物体1,2の運動量の合計としては

\[p=m_1v_1-m_2v_2\]

と表せます。これを使えば、物体1,2がぶつかってから物体2が壁にぶつかるという1サイクルでの運動量の変化を漸化式で表すことができるようになります。わかりやすく、サイクル開始時点での運動量を$p_n$、終了時点での運動量を$p_{n+1}$としてみましょう。

\[\begin{eqnarray}p_{n+1}&=&m_1\dfrac{p_n+m_2\sqrt{2k(\frac{1}{m_1}+\frac{1}{m_2})-\frac{1}{m_1m_2}p_n^2}}{m_1+m_2}-m_2\dfrac{p_n-m_1\sqrt{2k(\frac{1}{m_1}+\frac{1}{m_2})-\frac{1}{m_1m_2}p_n^2}}{m_1+m_2}\nonumber\\&=&\dfrac{(m_1-m_2)p_n+2m_1m_2\sqrt{2k(\frac{1}{m_1}+\frac{1}{m_2})-\frac{1}{m_1m_2}p_n^2}}{m_1+m_2}\nonumber\end{eqnarray}\]

さて、運動量$p$に関する漸化式ができました。この漸化式のnを数えれば衝突回数がわかるはずです(※1)。終了条件はもう衝突が起こらなくなるということで、それはすなわち運動量の絶対値が初期運動量の絶対値以上になったときと言えます(※2)。

※1：nが1増えるたびに衝突は2回起こります。
※2：これは正確には壁にぶつからなくなる条件で、この後に物体同士でぶつかるかは検証する必要があります。

シミュレーション

ここまで来たので、今までの計算が正しいことを確認するために実際にプログラムを書いて漸化式を数値的に解いてみます。

double n = 7;
double v1 = -1;
double v2 = 0;
double m1 = Math.Pow(10, n * 2);
double m2 = 1;

double p = m1 * v1 + m2 * v2;
double p0abs = Math.Abs(p);
double k = (m1 * v1 * v1 + m2 * v2 * v2) / 2;

uint collision = 0;

while(true) {
	var pprev = p;
	var sqrt = CalcSqrt(p, k, m1, m2);
	p = CalcNextMomentum(sqrt, p, m1, m2);
	collision += 2;

	if(Math.Abs(p) >= p0abs) {  // 合計運動量が初期運動量の反対向きを超えていたらもう壁にはぶつからない
		// -v2 > v1 ⇔ v1 + v2 < 0 だったらもう1回物体同士でぶつかる
		if(2 * pprev + (m2 - m1) * sqrt < 0)
			collision++;
		break;
	}
}

Console.WriteLine($"Pi = {collision / Math.Pow(10, n)}");


static double CalcSqrt(double p, double k, double m1, double m2)
{
	return Math.Sqrt(2 * k * (1 / m1 + 1 / m2) - 1 / (m1 * m2) * p * p);
}

static double CalcNextMomentum(double sqrt, double p, double m1, double m2)
{
	return ((m1 - m2) * p + 2 * m1 * m2 * sqrt) / (m1 + m2);
}

余談ですがC#10.0/.NET6.0になってnamespaceやMainメソッドが省略されるようになりました。このままのコピペでコンパイルが通ります。

運動量の漸化式を計算するメソッドがCalcNextMomentumですが、平方根部分のみを別途計算するCalcSqrtメソッドを用意しています。運動量を計算していって運動量の絶対値が初期運動量の絶対値以上となったらもうこれ以上の壁への衝突は発生しないものとして、物体同士で最後にぶつかるかどうかの判定へ移ります。

この時点での2つの物体の速度は、1回前のループで計算した運動量に対して計算した$v_1$と$-v_2$になるはずです。その大小関係ですので分子のみの差をとって判定しています。ここでルート部分が再利用できるため、メソッドを分けてわざわざルート部分の値を保存していたのです。

このプログラムは小数第7位までの計算（$n=7$）まで正常に動きます。それ以上は浮動小数点型の有効桁数が足りずに正確に計算できず、いつまでたってもループを抜け出せなくなってしまいます。

漸化式を解く

この漸化式解けるのか…？

まとめ

というわけで、多分私がこの問題だけを見て答えを見ていなかったら、ここまで解いて詰んで終わっていたでしょう。あの「エネルギー保存則を円の方程式、運動量保存則を直線の方程式とみたてて平面上に表すと、衝突回数が円周を運動量保存則の直線の傾き角の2倍で割ったものと考えることができる」とかいう天才的なひらめき、いったいどういう練習をしたら思いつけるようになるんだか。