由铁块碰撞出来的圆周率计算方法

作者：亚马逊的蝴蝶（Butterfly_of_Amazon）

---

果壳上有一篇文章《这一切要从碰撞的滑块和量子搜索讲起》，指出圆周率π与滑块碰撞次数有种奇妙的关系。今天给大家说说为什么会有这样的奇妙关系。

一、奇妙的关系

在光滑的平面上放着两个1 kg的滑块，在它们的左侧有一面坚固且不可移动的墙。两个滑块一个离墙近一些，一个远一些，而且后者正在向前者滑去。可以想象，接下来将会产生许多次碰撞。假设每次碰撞时都没有动能损失，整个过程将像下面图中展示的那样：

图1

右侧的滑块碰撞左侧的滑块，并将全部动量转移给左侧滑块，左侧滑块碰到墙反弹回来，再与右侧滑块发生碰撞并将动量全部转移过去。期间左侧滑块一共发生了3次碰撞。如果我们增大右侧滑块的质量，会发现有趣的事情：将右侧滑块质量换成100 kg，那每次碰撞时只会有一小部分动量发生转移，将会碰撞31次；如果将右侧滑块质量换成10000kg，那整个过程将发生314次碰撞。

不断100倍100倍地增加右侧滑块质量，你会发现神奇的现象：总碰撞次数与π每一位的数字越来越接近。

图2

二、这是真的吗

是不是很神奇？我通过计算，验证100kg、10000kg和1000000kg的情况，发现都与这个结论相符。验证过程如下：

图3

设滑块a和滑块b的质量分别为m_a和m_b，速度为v_a和v_b。速度随着碰撞而改变，设两个滑块间发生第n次碰撞后的速度分别为v_a,_n和v_b,n。

由于碰撞过程中没有动能损失，因此滑块a碰撞左侧墙壁后，速度将大小不变，方向相反。结合动量守恒定律，可列下式：

-m_av_a,n-1 + m_bv_b,n-1 = m_av_a,n + m_bv_b,n ........①

因为碰撞过程中总动能不变，可列下式：

m_av²_a,n-1 + m_bv²_b,n-1 = m_av²_a,n + m_bv²_b,n ...②

设m_b = K·m_a（图2中K为1000000，图3中K为100），由式①、②计算可得：

式③ 式④

当v_b,n为负，且|v_a,n|<=|v_b,n|时，两个滑块将不再发生碰撞，此时如果v_a,n>0，则滑块a还将与墙壁碰撞一次。将式③、④写成Excel表格中公式，设v_b初始速度为1，通过拖拽可方便地得到下表（表中以K=100为例）：

以表1中v_b为纵坐标，v_a为横坐标，可得下面半圆图形：

图4

细心的你可能已经看出：上图并不是真正的半圆，因为图中横坐标最大值是10，而纵坐标最大值是1，需要把纵坐标乘以10才是真正的半圆。为了方便描述，姑且把它称为半圆。

三、为什么会这样

为什么圆周率π与滑块碰撞次数之间有这样奇妙的关系？计算告诉你答案。

先问一个问题：根据表1中的数据，如何估算图4半圆的上半部分（也就是四分之一个圆）的面积S？

可能你已经想到了：可以把这四分之一个圆近似看成由8（约为n的最大值的一半）个竖长条组成，每个竖条的宽是相邻两点间的横坐标距离，也就是相邻两次碰撞后滑块a的速度差，高是这两点纵坐标的平均值，也就是相邻两次碰撞后滑块b的平均速度。将这8个竖条的面积之和当作S的近似值，虽然与真实值存在偏差，但通过增加圆上点的数量，也就是增加碰撞数量，从而增加竖条数量，可以减小偏差。

设碰撞次数的最大值为N，可列下式：

式⑤

再换个角度：把S近似为8个横条的面积之和，每个横条的高是相邻两点间的纵坐标距离，宽是这两点横坐标的平均值。可列下式：

式⑥

⑤+⑥，得：

式⑦

将③、④代入⑦，得：

式⑧

因为式②、m_b=K×m_a、v_b初始值=1，所以

式⑨ 故，式⑧ => 式⑩

由式⑨可知，所有坐标为 (v_a,√K·v_b) 的点位于一个半径为 √K 的圆上，此圆的面积为 K×π ；而S为此圆纵坐标被压扁到 1/√K 后面积的四分之一，因此有：

将⑩代入S，并考虑到在K非常大时，K/(K+1)≈1，可得：

N是两个滑块在前半程中相互碰撞的次数，后半程两个滑块相互碰撞次数也为N，滑块a与墙壁碰撞的次数为2N，故4N的整数部分即为总碰撞次数。这就是为什么圆周率π与滑块碰撞次数之间有这样的奇妙关系！

前面出现的式⑤提供了一个迭代计算圆周率π的方法。我通过迭代计算785次矩形面积，将圆周率精确到了小数点后5位（3.14159），计算过程中借助了Excel表格，有兴趣的朋友可以试试。