关于“用理想弹性碰撞能用来计算π”视频的小讨论
最近我在很多地方都看到了上面这张图，于是想拿出来分(蹭)享(个)一(热)下(度)，分享之前先摆上油管UP主3Blue1Brown的2个视频，其在16号以视频形式在油管上提出了这个问题并在20号公布了答案。其实早在1995年，美国东伊利诺大学的数学家Gregory Galperin就已经提出了这个问题，并在2003年为此写了一篇论文(上面的视频末尾提到了这个论文里提到的两面反光镜反射光线的内容)。后来在2014年，斯坦福大学的Gary Antonick提出的另一种比较直观的偏向于几何的证明方法，同时这也是上面视频里讲到的证明思路，这里呢就简单的引用一下2016年由“安安以迁迁”翻译的关于这两篇论文的内容，来简单说明一下。
想象地面上有一堵墙，墙的右边某处有一物体A，它的右边某处又有一物体B。假设地面无限长无限光滑，AB两物体都可视为质点，AB之间，以及它们和墙之间都作完全弹性碰撞。
先假设A和B的质量相同，A和B一开始都静止。我们朝左推一下B，让它有个初始速度。B运动一段时间后撞到A（第一次），由于动能守恒和动量守恒的缘故，B完全停止，而A以原先B的速度向左运动。A撞到墙（第二次）后以原速率向右弹回，运动后又撞到B（第三次）。接下去A完全停止，B向右一直运动。我们看到这样一共会发生3次碰撞。
为什么是3次？如果说这是因为圆周率π的第一个数字是3，你会不会觉得这种关联太牵强？让我们继续看。
保持其他假设不变，只是现在B的质量是A的100倍。最初B运动一段时间后撞到A，这时因为B的质量大于A，B并不会完全停止，而是会继续向左以比原先稍慢的速率运动，而A则会以比B更快的速率向左运动，直到撞上墙向右弹回，然后又撞上B向左弹回。这时B向左运动的速率又慢了一点，而A的速率则变得更快。这样A在B和墙之间反复碰撞，直至把B推动至向右运动。然后A仍在B和墙之间碰撞，B向右的速率将逐次加快，而A的速率则会逐次减慢，直至最终A赶不上B而不再发生碰撞为止。在下图的模拟中我们看到，这样一共会发生31次碰撞。
为什么是31次？如果说这是因为圆周率π去掉小数点后的前两个数字是3和1，你会不会觉得这纯属巧合？让我们继续看。
仍保持其他假设不变，只是现在B的质量是A的10000倍。下图模拟中的碰撞已变得看不清，但是通过计算，我们知道，一共会发生314次碰撞。
而假设B的质量是A的1000000倍，则会发生3141次碰撞；假设B的质量是A的100000000倍，则会发生31415次碰撞……总而言之，如果B的质量是A的10²ⁿ倍，那么会发生的碰撞次数就是把圆周率π的小数点拿掉后前面n+1位数表示的数字。只要两个物体一堵墙，推上一下，我们就能计算圆周率π。这显然不能再用巧合来解释了。
当然，在现实情况下想用这种办法稍微精确地计算一下π也是不可能的。首先，“无限光滑”“完全弹性碰撞”的条件太理想化。而A和B的质量不可能相差太大：比如n=50时，B和A的质量比要达到10¹ºº，即便A的质量小如电子，B的质量也将超过目前可观察宇宙的总质量。就不要说在极其微小和巨大质量的条件下，必须考虑诸如分子热运动，万有引力，量子效应，相对论效应等等物理学家们俯拾皆是的抬杠理由。我们这里谈论的与其说是一个物理问题，毋宁说是一个数学问题，也就是仅考虑经典的牛顿三大运动定律图景下的数学结论。
1995年美国东伊利诺大学的数学家Gregory Galperin在做一次关于小球碰撞的数学报告前发现了上面这个有趣的结论。当他在报告中公布这个发现时，开始听众们都觉得难以置信。但在给出证明后，听众们又纷纷表示信服。后来Galperin在2003年为此写了一篇论文。后面我要介绍的，则是由斯坦福大学的Gary Antonick提出的另一种偏向于几何的证明方法，比较直观。
先来看物理部分。
假设物体A质量为m，物体B质量为M，令k为M/m的（正）平方根，也即k²=M/m。为了一般性起见，我们不要求k是10的幂，甚至不要求k是整数，只要求k≥1，也即A的质量不大于B。在整个过程中，物体A的速度v(t)和物体B的速度V(t)都是时间t的函数，我们令向左的速度为正速度，这样向右的速度就被表示为负的。假设物体B初始的（向左的）速度为W。整个过程中动能守恒，于是在任意时刻t我们都有
1/2MV(t)²+1/2mv(t)² = 1/2MW²
也即
(V(t)/W)²+(v(t)/(kW))² = 1
这意味着在任何时刻t，点(V(t), v(t))都处于一个固定的椭圆上。如果我们作一个拉伸变换，定义新的关于t的函数x(t)=V(t)/W，y(t)=v(t)/(kW)，上面的公式就变成了
x(t)²+y(t)² = 1
也就是说不同时刻t所对应的p(t)=(x(t), y(t))都处在以原点为中心，以1为半径的圆上。
我们注意到，当在某段时间t₁到t₂之间没有发生过碰撞，那么这段时间里的v(t)和V(t)都是不变的，于是这段时间里点p(t)也是不变的。在整个过程中只会发生有限次碰撞，所以p(t)在坐标系中的图像将只有有限个点。下面我们将考察在某次碰撞前后，点p(t)会怎样变化。
如果是物体A和墙的碰撞，那么在碰撞前A向左运动，碰撞后则向右运动，而速度大小相同；物体B则保持原来的运动状态。这意味着在碰撞前p(t)处于x轴上方，碰撞后则处于x轴下方，和原来的点以x轴对称。
如果是物体A和B的碰撞，那么在碰撞前A必定静止或是在向右运动，不存在A和B都向左运动且B的速度大于A而与A相撞的可能（证明略，只指出一点：如果A正在向左运动，表明前一次碰撞发生在A和B之间），这意味着在碰撞前p(t)处于x轴上或是在它的下方。碰撞前后动量守恒，即
MV(t)+mv(t) = C
其中C在碰撞前后是个常数。如果将上式两边都除以MW/k，我们得到
kV(t)/W+v(t)/(kW) = kC/(MW)
也即碰撞前后
kx(t)+y(t) = kC/(MW)
上式等式的右边还是一个常数。这意味着碰撞前后的p(t)在同一条斜率为-k的直线上。综合上述结论，碰撞前后的p(t)的图像变化应该如下图所示，其中虚线的斜率为-k。
于是我们就得到了整个过程中p(t)变换模式：从点(1, 0)开始（对应着B刚开始以速度W向左运动，还未撞上A时的状态），交替地以斜率为-k的直线和以平行于y轴的直线成之字形在圆周上截出的点。下图是k=3时的情形。每一条蓝色虚线段代表一次碰撞，箭头方向则代表了此次碰撞前后系统状态（即点p(t)）转换的方向。
于是碰撞次数就是蓝色虚线段的数量，或者说是圆周上的p(t)图像的点数减1。这样，一个物理问题被转化为一个几何问题：从点(1, 0)开始，按照上述方法，交替地以斜率为-k的直线和以平行于y轴的直线成之字形在圆周上截出的点会有几个？而这个问题，即便从直觉上也可以感到它会和圆周率π很有关系。
顺便要指出的是，我们发现转换后的问题仅和k的值有关，或者说仅和物体B和物体A的质量比有关，而和物体A和物体B的具体质量、初始的推动速度W或是物体A和墙的距离等值无关。
在上节中我们把原来的物理问题转化成了一个几何问题：从点(1, 0)开始，交替地以斜率为-k的直线和以平行于y轴的直线成之字形在半径为1的圆周上截出的点会有几个？本节则将用几何的方法来解决此问题。
首先我们要解决的问题是，按照上述交替截圆周的方式，最后会以什么样的方式结束？也就是说，最后一个点的位置会是什么样的？答案是，这个点是最后一点当且仅当它落在下面这条红色的弧上，包括两个端点。其中红色虚线的斜率为1/k；蓝线为圆周在Q点的切线，所以它垂直于红色虚线，斜率为-k。因为圆半径为1，红弧的长度等于它所对应的圆心角大小，为arctg(1/k)。
这一点可以从物理和几何两方面来论证。
从物理图景来解释的话，不再发生碰撞的充要条件是物体A和物体B均处于静止或向右运动的状态，而且A的速率小于等于B的速率，即V(t)≤0，v(t)≤0且V(t)≤v(t)。翻译成函数x(t)和y(t)的语言，就是x(t) ≤ ky(t) ≤ 0，对应到圆周上，就是上图里的红弧。
从几何上说，如果圆周上的一个点p处于x轴下方或等于点(1, 0)，而且它不在上面那条红弧上，那么通过它的斜率为-k的直线必定处于上图蓝色切线的右上方，会在圆周上截出一个处于p的左上方的点p'来。换句话说，p必定不会是之字形截点过程的终结点。更进一步，如果上面这个新点p'不在红弧上（此时它会是截点过程的终结点），那么它必处于x轴的上方。同样，如果圆周上的一个点p在x轴上方，那么通过它的和y轴平行的直线必定可以在圆周上截出一个在x轴下方的点来，p也不会是之字形截点过程的终结点。
通过归纳法，从(1,0)出发的之字形截点过程必然可以一直持续下去，直到到达上图红弧上的一点。而一旦到达红弧，就既不可能用斜率为-k的直线向左上方截点，也不可能用和y轴平行的直线向下截点：它必然是终结点。
我们看到物理和几何的方法都得到了同样的结论：之字形截点过程终结，当且仅当截点出现在上述的红弧上。
而截出终结点的方式有两种，下面我们分别讨论。
第一种方式是终结点被斜率为-k的直线所截出。在物理上对应的过程是倒数第二次物体A和墙相撞向右弹回，并追上正向右运动的物体B，进行了最后一次撞击；撞击后物体A和B仍朝右运动，只是A的速率小于等于B。第一节中质量比为100的动画中的撞击过程就是这样结束的。在k=3时终结点也是这样被截出，下面我们以此为例说明。
在图像中我们将每个点标了号，并通过初始点1作了切线。这样除了最后一点10外，每一点都对应着以它为顶点，由两条蓝色射线形成的角。因为所有斜率为-k的直线都平行，所有平行于y轴的直线都平行，于是角1到角9依次两两都是平行直线的内错角，它们均相等，其正切是1/k，故等于arctg(1/k)。它们是圆周角（角1是弦切角，可以看作是圆周角的特殊情况）。从中学平面几何里我们知道，它们对应的圆弧长度均相等；因为圆半径为1，它们所对应的弧长是两倍的圆周角大小，即2arctg(1/k)。在k=3时，弧1-2、弧1-3、弧2-4、弧3-5、弧4-6、弧5-7、弧6-8、弧7-9，弧8-10的长度均是2arctg(1/3)；而弧9-10的长度小于此值。出现最后这段弧长会等于前面的弧长的情况仅当终结点正好是点(-1, 0)，这是k=1也就是物体A和物体B质量相等的情况。
第二种方式是终结点被平行于y轴的直线所截出。在物理上对应的过程是倒数第二次物体A与物体B相撞后向左运动，最后撞在墙上向右弹回，但此时它的速率已小于B的速率，再也追不上B。这是第一节中质量比为10000的动画中的撞击过程结束的方式，同样也是k=4时的情况。下面我们以k=4为例说明。
和k=3的情况很相似，圆上的截点把圆周截成了若干段长度为2arctg(1/k)的弧，以及最后一段弧12-13。最后这段弧的长度总是小于等于两倍的红弧（因为最后一点须在红弧内），也就是2arctg(1/k)。
综合上述两种方式，我们就得到结论：之字形在半径为1的圆周上截出的点的个数，正是圆周的长度除以2arctg(1/k)的结果的整数部分再加1，这个1对应着那段比较短的弧。（严格地说，如果圆周的长度除以2arctg(1/k)的结果恰好是整数，如当k=1时，那么就不用加这个1了，这一点将在后面补充说明）。或者说，碰撞的次数（它等于点的个数减1）就是圆周的长度除以2arctg(1/k)的结果的整数部分。但半径为1的圆周的长度正是2π。这一下，碰撞次数和圆周率的关系可谓昭然若揭，剩下的只不过是一点计算细节罢了。
碰撞的次数为圆周长2π除以2arctan(1/k)的整数部分，也就是[π/arctan(1/k)]，其中[]代表取整运算。当k大于等于1时，arctan(1/k)总比1/k大一点，但差不了多少；更精确地，π/arctan(1/k)和π/(1/k)=kπ差不多大小，k越大，这个差就越小。这一点可对f(x)=π/arctan(x)-π/x在零点附近作泰勒展开并作估计即可，这要用到数学分析知识，这里就不具体进行了。于是[π/arctan(1/k)]和[π/(1/k)]=[kπ]一般来说是相等的。
总的来说，在理论上[π/arctan(1/k)]和[kπ]有极小的可能差1。一种可能是上述的圆周的长度除以2arctg(1/k)的结果恰好是整数的情况，那么很好，因为arctan(1/k)总比1/k大一点，于是[kπ]等于[π/arctan(1/k)]-1，正是圆周上的点个数减1。
另一种可能则相当麻烦。如果π/arctan(1/k)和kπ恰好在某个整数的两侧，就如2.00001和1.99999在2的两侧那样，相差极少，取整结果却差1。为了对这个细节做精确估计，Galperin差不多用了论文的三分之一篇幅来作相当繁琐的讨论，也没有完全排除这种可能性。但对于我们来说，即便[π/arctan(1/k)]和[kπ]真的差1也不是太有所谓，因为这无非是在说，碰撞的次数和π的前几位数几乎一样，最多差1。这并不影响结果的漂亮，也不影响我们对为何初始的碰撞问题会和圆周率有出人意料却在数理之中的联系的理解。所以我们接下来就当这种可能性不存在。
这样我们就知道，总的碰撞次数是[kπ]。于是最终结论就显而易见了，在k=10时碰撞次数是[10π]，k=100时碰撞次数是[100π]，等等等等，就是圆周率π的前几个数字。可以看出，如果不用十进制而是用其他进制，也会有类似的结果。比如说在二进制下，圆周率π约等于11.0010010000111111……，如果我们取k=2⁵=32，也就是质量比为2¹º=1024时，模拟结果表明会发生100次碰撞：
而100在二进制下的表示为1100100，恰是π的二进制表示的前7个数字。
参考文献：
[1] Gary Antonick, The Pi Machine, Numberplay, Wordplay Blog, http://wordplay.blogs.nytimes.com/2014/03/10/pi/
[2] Gregory Galperin, Playing pool with π (the number π from a billiard point of view), Regular and Chaotic Dynamics, Volume 8, Number 4, Pages 375–394 (2003).
[3] 安安以迁迁, 碰撞出来的圆周率, 简书, https://www.jianshu.com/p/3dea48c0fc64, https://www.jianshu.com/p/b9d3d4918ca7, https://www.jianshu.com/p/bb5eb822db78
