算圆周率的奇怪方式增加了——Python利用物理定律模拟计算圆周率

用动量和能量守恒计算圆周率

实现原理

灵感来自G. Galperin的论文PLAYING POOL WITH π (THE NUMBER π FROM A BILLIARD POINT OF VIEW)，DOI: 10.1070/RD2003v008n04ABEH000252。我在这里只是简单介绍下原理，感兴趣的同学们可以自己仔细阅读。论文是全英文，有需要可以私信我，有问题我帮你解答。Youtube上3Blue1Brown的视频有介绍这片论文，有很好的图片和讲解，更好理解。

简单概括下，Galperin的想法很神奇。不同于普通的算π的方法（如无限级数法），作者用动量守恒和动能守恒算$\pi$！

首先，我们来理解一下动量守恒和能量守恒。
在完全弹性碰撞中，两个物体在碰撞后除了动能的转换没有其它能量的消耗。这是完全弹性碰撞的动量守恒公式：
$$m_1{v}_1+m_2{v}_2=m_1{v}_1^{{}^{\prime }}+m_2{v}_2^{{}^{\prime }}$$

基于此，二者能量的关系可用动能守恒表示：
$$\frac{1}{2}{m}_1{v}_1^2+\frac{1}{2}{m}_2{v}_2^2=\frac{1}{2}{m}_1{v}_1^{{}^{\prime }2}+\frac{1}{2}{m}_2{v}_2^{{}^{\prime }2}$$
连立上面两条公式，我们得到：
$$\{\begin{array}{lr}{m}_1{v}_1+m_2{v}_2=m_1{v}_1^{{}^{\prime }}+m_2{v}_2^{{}^{\prime }}& \\ & \\ \frac{1}{2}{m}_1{v}_1^2+\frac{1}{2}{m}_2{v}_2^2=\frac{1}{2}{m}_1{v}_1^{{}^{\prime }2}+\frac{1}{2}{m}_2{v}_2^{{}^{\prime }2}& \end{array}$$
化简该公式，我们得到：
$$\{\begin{array}{lr}{v}_1^{{}^{\prime }}=\frac{(m_1-m_2)v_1+2m_2{v}_2}{m_1+m_2}& \\ & \\ v_2^{{}^{\prime }}=\frac{(m_2-m_1)v_2+2m_1{v}_1}{m_1+m_2}& \end{array}$$
这样，我们就有了程序模拟的核心算法。

在下图所示的理想的物理系统（无摩擦力，墙和地面质量无限，球体无角动量等）中，这两条关系成立：（图源：Galperin）

具体的数学证明，Galperin在论文中有提到。我们要做的，就是用暴力模拟验证他的结论。设物体朝墙运动时速度为正，远离墙时速度为负。二者的碰撞在以下情况时：$$0>v_2\ge v_1$$，两个物体不会再次相遇。

Galperin的数学论证非常详尽。总而言之，他证明了以下关系：
$$\Pi (N)=314159265358979...$$
其中$\Pi (N)$是当大物体重量是小物体的$100^N$倍时的碰撞次数，$N$是计算$\pi$的位数-1。
我们要做的，就是根据他的论文写个程序实现。

代码

#!/usr/local/bin/python3
import time
from decimal import *
from decimal import Decimal as dcm #这个库提高下浮点运算精度

getcontext().prec = 50 #小数位数设为50位

N = int(input("Enter N:"))
m = dcm(1)        #初始化m1
M = dcm(m*100**N) #初始化m2

vm0 = dcm(0)      #初始化v1
vM0 = dcm(10000)  #初始化v2

vm = vm0 	#初始化v1'
vM = vM0	#初始化v2'

cnt = -3 #初始化计数器 经过测试软件总会多算几次 所以先-3

start = time.time() #设置一个计时器，算所需时间
while not (vm < vM and vM < dcm(0)): #也就是当两个物体不会同时向右且不再相遇时
    vm = dcm((m-M)*vm0+2*M*vM0)/dcm(M+m)
    vM = dcm((M-m)*vM0+2*m*vm0)/dcm(M+m) #上面提到的两个核心模拟算法
    cnt += 1
    vm0 = dcm(-vm) #弹性碰撞，小物体的动量变成反方向
    vM0 = dcm(vM)  #大物体动量更新，方向不变
    cnt += 1 #又撞了一次
end = time.time() #结束计时器
print("Number of collisions:", cnt)
print("Time elapse:", end-start, "seconds.") #输出

运行结果

运行结果见下表：（本人自制）
最左第一列是N的值
第二列是输出的碰撞次数
第三列是运行时间
第四列是根据输出计算的$\pi$值
第五列是和真实$\pi$值的差距

下图是N的变化与计算偏差的指数回归：（作者自制）

由上图可见，随着N的增大，计算的圆周率指数级地无限趋于准确。

下图是N的变化与计算所需时间的指数回归：（作者自制）

可见，随着准确度的指数级提高，计算的时间复杂度也指数级地提高。因此，一般的笔记本真算不过来太大的N。

总结与拓展

我们的程序确实非常有效。然而，有很多其它计算$\pi$的方法更加高效。我们的程序也有一些疏漏。比如，输出的值永远是奇数。这是因为在while循环里默认最小增加值是2。有兴趣的小伙伴可以思考下如何解决！也许可以多加几个判定？
对于Galperin论文感兴趣的同学请私信我或留评论，我会很乐意帮助你们解答相关问题！

ps：本人在用这些东西写小论文，所以转载请引用！