从碰撞到圆周率

· · 个人记录

Hihi大家好，这里是许久没有发博客的computer！

最近在学校学习到了牛顿力学最后一部分——动量的相关知识，忽然想起3Blue1Brown

大佬曾经讲述过如何从碰撞中获得\pi,那我今天也来炒一炒冷饭，来讲述如何简单地

写一个程序来计算\pi的前几位数！

事不宜迟，我们开始吧！

前置知识

机械能守恒定律 在只有重力和弹力做功的系统中，系统机械能守恒。

动量守恒定律 如果一个系统不受外力或者所受外力矢量和为0，那么系统动量守恒。

讲解

1、证明结论

在这里先贴上3Blue1Brown大佬的原视频链接：

link

也可以看我的表述：

如图所示，在光滑水平面上有m_1,m_2两个物块，左侧有一堵墙，m_2有向左的初速度v_2。

假设全部过程为弹性碰撞，m_1撞击墙壁后速度只改变方向。

请证明：当\sqrt{\tfrac{m_1}{m_2}} \rightarrow 0时，m_1,m_2的撞击次数n \rightarrow \pi \sqrt{\frac{m_2}{m_1}}.

即:求证\lim\limits_{\sqrt{\tfrac{m_1}{m_2}} \to 0} n = \pi \sqrt{\frac{m_2}{m_1}}.

证明： 设水平向右为正方向。

对于m_1,m_2组成的系统，由机械能守恒和动量守恒可得

&p = m_1v_1 + m_2v_2 & \\ &E_k = \tfrac{1}{2}m_1v_{1}^{2} + \tfrac{1}{2}m_2v_{2}^{2}& \\ \end{matrix}\right.

并且p,E_k为定值。

很明显这两个方程很有解析几何的味道，但是别急，要先对式子进行一些变换。

设x=\sqrt{m_1}v_1,y=\sqrt{m_2}v_2,则方程转化为

&p=\sqrt{m_1}x + \sqrt{m_2}y & \\ &x^2 + y^2 = 2E_k & \\ \end{matrix}\right.

这时候很明显动量的式子变成了一条直线的解析式，动能的式子变成了一个圆。

同时动量解析式的斜率k=-\sqrt{\tfrac{m_1}{m_2}},tan\theta =\sqrt{\tfrac{m_1}{m_2}}.

此时我们再来把方程放进平面直角坐标系中，可以看到：

首先是一个半径为\sqrt{2E_k}的圆，表示为动能；一条直线经过动能圆与y轴负半轴交点，

表示为动量，其中斜率满足上述条件。

由于初始时v_1=0,v_2=v<0，所以自然直线相交于动能圆与y轴负半轴交点。

因为整个过程中的碰撞满足动量守恒和机械能守恒方程，所以圆与直线的交点就表示一个碰撞状态。

接下来m_1撞击墙壁导致速度反向，所以将会使直线平移到y轴对称点位置。

同理，接下来m_1速度方向改变后将会出现一条平行于第一条斜线的直线：

同理，我们可以一直画到不能再画为止：

很明显这个时候再画的话就与圆没有交点了。

由于所有斜线平行，所有的水平线平行，所以我们可以标一下角度\theta：

可以看到所有的\theta都对着一段圆弧，而我们的目标是要用尽量多的圆弧覆盖整个圆周。

又因为\theta为这些圆弧的圆周角，同弧所对圆周角等于圆心角一半，

那么本质上来说就是说这个问题转变为了

求多少个2\theta角所对圆弧能够覆盖整个圆弧。

可以稍微停下来思考消化一下

即为：

求多少个\theta角所对圆弧能够覆盖整个圆弧一半。

那么就有如下方程：

&n\theta < \pi & \\ &(n+1)\theta > \pi & \\ \end{matrix}\right.,n\in \mathbb{R}

解之得：\tfrac{\pi}{\theta}-1\leqslant n < \tfrac{\pi}{\theta}.

当\sqrt{\tfrac{m_1}{m_2}} \rightarrow 0时

\tan\theta\approx\theta\approx \sqrt{\tfrac{m_1}{m_2}}

故可知：\lim\limits_{\sqrt{\tfrac{m_1}{m_2}} \to 0} n = \pi \sqrt{\frac{m_2}{m_1}}.

证毕。

2、计算\pi

实际上计算碰撞次数非常简单。

首先我们回到这个方程，并进行一定的修改：

设v_1,v_2为m_1,m_2碰撞前的速度，v_1\prime,v_2\prime为m_1,m_2碰撞后的速度，

故可知：

&m_1v_1 + m_2v_2 = m_1v_1\prime + m_2v_2\prime& \\ &\tfrac{1}{2}m_1v_{1}^{2} + \tfrac{1}{2}m_2v_{2}^{2} = \tfrac{1}{2}m_1v_{1\prime}^{2} + \tfrac{1}{2}m_2v_{2\prime}^{2}& \\ \end{matrix}\right.

通过一顿细心而又复杂的计算，我们可以得出以下结论：

&v_1\prime = \tfrac{(m_1-m_2)v_1+2m_2v_2}{m_1+m_2} & \\ &v_2\prime = \tfrac{(m_2-m_1)v_2+2m_1v_1}{m_1+m_2} & \\ \end{matrix}\right.

这个计算应该是正常操作了吧。

再根据分析可知，如果v_1<0，那么m_1一定会再和墙碰撞一次；

如果v_1>v_2（思考一下，需要考虑方向吗？），那个m_1会与m_2碰撞一次，

然后速度满足上述结论。

只有当0\leqslant v_1 < v_2时不会增加碰撞次数，结束计算。

综上所述，我们可以写出如下代码：

#include <bits/stdc++.h>
using namespace std;
long double m1 = 1,m2,v1 = 0,v2 = -10,v_1,v_2;
int main()
{
    long double rate,cnt = 0;
    scanf("%Lf",&rate);
    long double div = rate;
    rate *= rate;
    m2 = m1 * rate;
    while(1)
    {
//      printf("%lf %lf\n",v1,v2);
        if(v1 >= 0 && v2 > v1)
            break;
        if(v1 < 0)
        {
            cnt++;
            v1 *= -1;
        }
        if(v1 > v2)
        {
            cnt++;
            int tmp = v1;
            v1 = ((m1 - m2) * v1 + 2.0 * m2 * v2) / (m1 + m2);
            v2 = ((m2 - m1) * v2 + 2.0 * m1 * tmp) / (m1 + m2); 
        }
    }
    printf("%Lf",cnt/div);
    return 0;
}

使用时输入10^n,n\in \mathbb{N}的数字来获得\pi的小数点后前n为数。

例如输入"1000000",获得"3.141592"。（需要注意不带四舍五入！）

需要注意的是，我用了C99才加入的long int变量，编译前要在Dev-C++的

“工具”——“编译选项”——“代码生成/优化”中将“支持所有ANSI C标准”选择“Yes”即可。

总结

这个方法我觉得还可以，可惜在计算8位的时候已经要花费6s左右，9位时候就计算了很久就没下文了。

但不管怎么说，通过物理模型计算\pi都是不同于以往数学的新颖方法，

这告诉我们自然科学的背后本质是数学，许多看似平常的物理问题背后是数学在发挥着作用。

参考资料：

（1）3Blue1Brown,So why do colliding blocks compute \pi?