木块碰撞次数暗藏圆周率：用Python模拟复现3.1415

木块碰撞和π的联系是我这几年见过最反直觉、也最值得亲手验证一遍的物理现象。光滑地面上放一大一小两个木块左边一堵刚性墙把小木块夹在大木块和墙之间给大木块一个向左的初速度剩下的就是数碰撞次数大木块撞小木块、小木块撞墙、小木块再撞大木块……当你精心调整两个木块的质量比统计出来的碰撞总次数会像变魔术一样逐位吐出圆周率π的数字。这个现象适合谁物理爱好者会喜欢它背后的相空间解释编程学习者可以把它当成一个绝佳的模拟练习物理老师能直接拿来做课堂演示参加科创项目的同学也能从中挖出做实验的思路。我自己第一次看到结果时以为程序写错了后来把原理啃透才发现离散的碰撞计数和连续的圆周率之间确实有一条严谨的几何桥梁。下面从现象、原理、程序模拟、实物实验到避坑经验完整讲一遍。1. 这个实验到底在验证什么1.1 一个一眼看上去像脑筋急转弯的场景先把场景交代清楚。地面绝对光滑墙绝对刚性两个木块和墙的碰撞都是完全弹性没有能量损失。大木块质量记为M小木块质量记为m小木块放在大木块和墙之间初始让大木块以某个速度向左运动小木块静止。所谓“一次碰撞”指的是任意两个物体之间接触并分离的一次事件包括大木块撞小木块、小木块撞墙、小木块再次撞大木块。如果M和m相等全过程一共3次碰撞。你可以自己演算一遍大木块撞上静止的小木块因为质量相等速度交换大木块停下小木块向左撞墙撞墙后小木块反弹向右再撞上静止的大木块再次交换速度小木块停下大木块向右离开。整个过程3次撞击不多不少。把质量比改大一点奇迹就出现了。M:m 100:1时碰撞总次数是31M:m 10000:1时碰撞总次数是314M:m 1000000:1时碰撞总次数是3141。我列出几个常用质量比和对应的碰撞次数质量比 M:m模拟出的碰撞次数看到的数字1:133100:1313.110000:13143.141000000:131413.141100000000:1314153.1415表格里的数字已经说明问题质量比每放大100倍碰撞次数的数字序列就相当于π的位数向后多显现一位。离散的碰撞次数竟然能编码一个无限不循环小数这放在任何领域都很反直觉。1.2 为什么碰撞次数会逐位吐出圆周率最先冒出来的念头肯定是“巧合”。但稍微验证就会发现不是。把质量比改成别的不规则数值比如500:1碰撞次数是一个没有明显规律的整数并不会刚好对应某个熟悉的常数。只有把质量比取成100的整数次幂也就是100、10000、1000000……的时候输出才会呈现出3、31、314、3141这样的序列。原因是质量比经过开方之后直接决定了碰撞在某个抽象空间里“每步转过的角度”。M:m 100的k次方对应√(m/M) 10的负k次方。因为这个比值极度接近0每一步转过的圆心角也极其接近10的负k次方弧度。而整个过程里状态点在那个抽象空间里大约走完半个圆周也就是转过π弧度。于是碰撞次数近似等于π除以每一步的小角度得到π乘以10的k次方取整后正好就是π从小数点后逐位展开的数字。简洁说质量比决定了“步长”墙和弹性碰撞决定了“总路程”π是那个总路程里天然存在的几何常数。木块并没有偷偷计算π它只是在一个圆形几何结构中走了很多很多步步数正好携带了π的信息。2. 看懂原理先看“速度空间”2.1 把两个速度变到一个圆上要深入理解这个现象最优雅的办法是不要盯着物理空间里的坐标看而是把两个木块的速度放到一个抽象的“速度空间”里。设大木块速度为vA小木块速度为vB。完全弹性碰撞过程中系统满足两条守恒律动量守恒和动能守恒。动能守恒的表达式是0.5M vA² 0.5m vB² E这个方程在普通的vA-vB坐标系里是一个椭圆不太好用。如果做一个坐标变换令横坐标P √M · vA纵坐标Q √m · vB那么动能守恒就变成P² Q² 2E。这是一个标准的圆。动量守恒则写成√M · P √m · Q 常数。在(P, Q)平面上这是一条直线。系统的任意一个状态也就是两个滑块任意一组速度组合都必须同时满足“在圆上”和“在直线上”。如果没有碰撞状态点是不动的一旦发生碰撞能量和动量仍然守恒但速度发生突变于是状态点沿着这条直线与圆的交点跳到另一个位置。2.2 碰撞就是圆圈上的转跳接下来看两类碰撞分别做了什么。小木块撞墙时它的速度vB直接反向也就是Q变为-Q。在速度空间里这相当于状态点关于横轴做一次镜像反射。大木块和小木块互相碰撞时动量守恒方程那条直线的方向由质量比决定碰撞后状态点会跳到这条直线与圆的另一个交点上本质上也是一次反射。两种反射叠加状态点就在圆周上离散地跳来跳去每一步跳过的圆心角只由质量比决定。当质量比极大时√(m/M)是一个非常小的数状态点每一步跳过的角度也极小。墙在左边意味着整个过程的状态点不会无限跳下去它有一个明确的终点大木块不再追上小木块小木块也不会再撞墙。从初始点跳到终点状态点在圆周上扫过的总角度趋近于π。这就是为什么碰撞次数大致等于π除以步长而步长又约等于√(m/M)。形象地说质量比越大状态点在圆上的“脚步”越小走完同样一段半圆弧需要的步数就越多。每100倍的质量比相当于步子缩小到原来的十分之一步数就多出十倍。再加上一个π的系数就得到了3、31、314、3141这样逐位展开的数字。3. 程序模拟从零开始造一台“π打印机”3.1 物理模型与事件驱动循环原理说得再漂亮不如让代码跑一遍来得踏实。自己写模拟器别看答案看它能不能真的输出314。物理模型很直接墙固定在x0处向右为正方向。大木块初始在右侧小木块在大木块左侧、更靠近墙。大木块初速度为负即向左运动小木块初速度为零。每次循环只需要计算两个候选事件小木块什么时候撞墙以及两个木块什么时候互相碰撞。谁先发生就先推进到那个时刻执行对应的碰撞计数器加一。这里有个新手最容易踩的坑不能简单用“如果小木块速度小于0就处理撞墙如果大木块速度小于小木块速度就处理两木块碰撞”这种固定顺序判断。因为在某些瞬间两件事都可能即将发生必须比较它们的发生时间先处理最早发生的那个。我见过很多版本的程序把顺序写反结果初始状态下两个木块明明马上要相撞程序却先去处理小木块撞墙速度更新顺序整个乱掉。3.2 完整代码与运行结果我用分数来跑因为普通浮点数在质量比很大的时候会带来误差。分数运算慢一点但验证到1000000:13141次或者100000000:131415次完全足够。from fractions import Fraction # 两个滑块的质量改这两个数就能切换“精度” M Fraction(10000, 1) m Fraction(1, 1) # 墙在x0处滑块坐标向右为正 xA Fraction(10, 1) # 大滑块初始在右侧 xB Fraction(9, 1) # 小滑块初始在大滑块左侧、靠近墙 vA Fraction(-1, 1) # 大滑块初始向左运动 vB Fraction(0, 1) # 小滑块初始静止 hits 0 def advance(t): global xA, xB xA vA * t xB vB * t def next_events(): t_wall None if vB 0: t_wall -xB / vB t_ab None if xA xB and vA vB: t_ab (xB - xA) / (vA - vB) return t_wall, t_ab while True: t_wall, t_ab next_events() if t_wall is None and t_ab is None: break if t_ab is not None and (t_wall is None or t_ab t_wall): advance(t_ab) vA_new ((M - m) * vA 2 * m * vB) / (M m) vB_new ((m - M) * vB 2 * M * vA) / (M m) vA, vB vA_new, vB_new hits 1 else: advance(t_wall) vB -vB hits 1 print(hits)把M换成100、m保持1输出31M换成10000输出314M换成1000000输出3141M换成100000000输出31415。看到这些数字从自己写的代码里蹦出来那种感觉和看别人截图完全不一样。值得说明的是坐标初始值其实不影响碰撞次数只要大木块在小木块右侧、小木块在墙右侧就够了。位置只影响每次碰撞发生的具体时刻不影响事件的顺序和次数。真正决定结果的是质量比以及所有碰撞必须保持完全弹性这个前提。3.3 让程序输出更多π位想看到更长的π位序列可以继续增大质量比。M:m 10^10时程序会输出314159到10^12时浮点数开始不稳定而分数版本虽然精确但循环次数已经达到三亿次量级跑起来很吃力。这里提醒一句不要盲目追求“用模拟器打印出几十位π”因为每一步都是Fraction分数运算几千次还好几亿次会让运行时间变得很难看。想验证更多位正确做法是从几何原理出发推导闭式公式或者用高精度浮点库优化循环。作为演示项目验证到314159已经很能说明问题。4. 参数选择与误差控制4.1 质量比为什么要取100的幂这是整个验证里最容易被忽略的前提。质量比取100、10000、1000000是因为它们开方后是10、100、1000正好让√(m/M)成为10的负整数次幂。这个“步长”是十进制下的整齐数字碰撞次数才会是3、31、314、3141。如果你随手设成M:m 200:1√(m/M)约等于0.0707算出来的碰撞次数是一个看起来毫无规律的整数和π的关系要经过额外换算才能发现。不是说其他比例就没有数学规律而是它们不会直接输出十进制下的π数字。想向别人演示这个现象从100:1开始一步一步乘100是最好的路径。实操顺序我建议这样先跑1:1确认输出3再跑100:1确认输出31然后10000:1确认输出314。每一步都验证通过再往上加。如果一上来就跑10^8一旦输出不对你很难判断是原理理解错了、代码写错了还是浮点精度问题。4.2 浮点数、分数和终止条件写模拟器时大多数人第一版会用float。质量比在10^8以内float基本够用我实测到10^10时还能输出314159再往上就可能出现偏差。原因是接近某个循环末端时碰撞时间差会变得极小两个事件几乎同时发生浮点误差会把事件顺序颠倒导致少计一次或多计一次。解决办法就是换成Fraction或Decimal。Fraction的优势是完全精确没有任何舍入误差缺点是慢。Decimal需要自己设置精度位数但实现起来也简单。更隐蔽的问题是终止条件。我见过有人用“大木块速度大于0就停止”来结束循环这在某些参数下能蒙对但不是通用判断。正确的做法是持续计算下一次事件时间直到“没有墙碰撞事件也没有两木块碰撞事件”为止。我的代码里用next_events()返回两个时间一个是None就跳过对应事件两个都是None就退出循环这个结构在三种情况下都不会卡死。5. 实体实验可行吗5.1 气垫导轨方案看完程序输出很多人会想那我把木块放在桌上碰一下能不能也数出31次想法很好实体实验的难点要从大到小说出来。第一是摩擦普通桌面滑块跑不了几步就停了第二是非弹性两个木块碰撞后会损失动能速度越小碰撞次数越难数第三是计数几十次碰撞全靠人眼数容易眼花。所以实体实验最现实的方案是用气垫导轨。气垫导轨可以把摩擦降到很低滑块基本悬浮。挡板一端固定一个刚性端面当“墙”。小滑块靠近挡板放置大滑块放远一点给大滑块一个较小的初速度。质量比做到100:1也就是小滑块1份质量、大滑块配重到100份理论上能观察到31次碰撞。这个量级在气垫导轨上是可行的再往上到10000:1需要314次碰撞实验时间会拉长摩擦和气源稳定性都会成为大问题。我自己试过的经验是先不要追求一次成功把气垫导轨调平至少花十分钟滑块在任意位置都能保持静止或极慢匀速才算调好。其次初速度要小速度太大碰撞时挡板会震速度分布不均匀也会影响后续计数。5.2 计数方式与实操心得计数是实体实验最劝退的部分。人眼数31次已经有点累数314次完全不现实。推荐两种方式。第一种光电门。把光电门放在撞墙点附近小滑块每次经过挡板触发一次信号配合计数器直接累计撞墙次数。但两滑块之间的碰撞可能不经过光电门需要再在碰撞区域多放一个门两路信号合并计数。第二种高速摄像。现在的手机慢动作模式就能胜任拍完之后逐帧数碰撞次数。好处是可以反复回放不会漏数坏处是后期处理费时间。拍的时候记得让两个滑块的颜色区别明显背景干净后期数帧会轻松很多。还有一个经验实体实验别指望验证到3141次。气垫导轨在100:1能做到31次已经很有说服力了再往上属于“原理演示”而不是“工程复现”。想要更长的π数字把舞台交给程序模拟就好。6. 常见问题与排查技巧实录6.1 三个高频翻车现场第一个翻车点是循环顺序写错。现象是程序一开始就不对或者输出结果不稳定。解决办法是把两件事的“发生时间”算出来再比较而不是固定一个if顺序。我在3.1里已经强调过这里再重复一遍因为这是网上很多版本代码的隐藏毛病。第二个翻车点是浮点精度导致末尾数字偏一位。现象是质量比到10^10以后有时输出314159有时输出314158或者比预期少一次。解决办法是把float换成Fraction精度问题立刻消失。如果用的是Decimal记得把prec设置到50位以上因为中间过程会出现极小的速度差默认精度不够用。第三个翻车点是计数口径不一致。同样一套物理过程“3次”还是“4次”可能只是因为你把初始大滑块第一次撞小滑块算进去了或者把最后停下来的瞬间误判成一次碰撞。我的代码里每执行一次碰撞事件就加1事件由时间推进驱动不会出现模糊口径。6.2 问题速查表现象可能原因解决办法程序输出不是整数规律序列质量比没有取100的幂改用1:1、100:1、10000:1起步程序一开始就乱碰撞事件顺序固定判断错了计算并比较两个事件的时刻大质量比时末尾数字偏一位float浮点误差用Fraction或Decimal能量不守恒碰撞公式写错核对vA_new、vB_new公式实体实验次数波动大气垫导轨没调平、初速度太大先调平导轨用小初速度多次重复计数器多计/漏计光电门触发阈值设置不当用高速摄像复盘确认触发逻辑表格之外还有一个隐蔽问题代码里用FractionM和m必须是Fraction类型否则整数和Fraction混乘时会有类型提示和运算速度上的麻烦。我习惯把所有初始量都包一层Fraction再参与后续运算。7. 一个物理现象的延伸思考7.1 圆周率藏在旋转对称里这个现象最迷人的地方在于它把离散的碰撞计数和连续的圆周率联系起来而背后的根源是“圆”。能量守恒在变换后的坐标系里天然是一个圆守恒定律在碰撞过程中迫使状态点在圆周上移动π作为圆的周长与直径之比自然就跑到了计数结果里。类似的事情在物理里反复出现。圆周率出现在简谐运动周期里出现在电场和引力场的一些积分解里出现在量子力学波函数的相位因子里。凡是涉及旋转对称或者圆形相空间的问题π都有可能出现哪怕表面上只是两个木块在一条直线上互相撞来撞去。理解这一点之后再看“碰撞次数为什么是314”就一点不玄了。它不是木块有灵性而是数学结构在底层起作用。做这类验证项目的价值不只是得到一个数字而是养成“把物理量换到合适坐标系里再看问题”的习惯。7.2 这个项目还能怎么玩验证完基本现象还能往三个方向扩展。第一个方向是可视化。把每次碰撞后的(P, Q)状态点画在相平面上你会看到一串点在一个圆弧上跳动每跳一次就是一次碰撞。这个图比干看数字直观得多强烈建议自己画一次。第二个方向是换个碰撞规则。把完全弹性改成恢复系数小于1的情况观察碰撞次数和π数字如何退化能帮助理解“能量守恒”在现象中的关键作用。第三个方向是求闭式公式。既然碰撞次数近似等于π除以步长那么精确公式应该是一个包含反正切函数的取整表达式。对这个感兴趣的读者可以在模拟结果的基础上反推公式再和模拟结果对比这本身就是一个很有意思的数学小课题。最后分享一个我自己的体会验证这个项目时最难的部分不是写程序而是说服自己“这不是巧合”。把速度空间画出来、把每一步碰撞的角度量化、把100的幂和十进制展开联系起来这套思考过程比最终输出314那几个数字珍贵得多。建议每一位读者都亲手跑一遍最简单的1:1情况盯着一组速度变化看三分钟再去跑100:1你就能直观感觉到每一步“转”得更密集自然理解为什么碰撞次数会变多。如果只记住一个技巧那就是永远从质量比1:1开始调试确认代码和物理都正确后再放心地乘100往上加。这样即使后面输出和预期不一致你也能快速定位是精度问题还是逻辑问题而不是面对一团乱麻。