用Python模拟乒乓反弹：从固定步长到事件驱动

我第一次认真琢磨“乒乓原理”倒不是因为打乒乓球而是同事调试一段数据流的时候说了句“这不就是个乒乓缓冲嘛”。后来才发现这个词在不同圈子里指的东西完全不是一回事搞硬件和图形处理的说的是数据在两块缓冲区之间来回倒换搞物理模拟的说的却是一粒小球落到水面后反弹、下落、再反弹的振荡过程。这篇博文要聊的是后面这个版本也就是很多动力学入门题里出现过的那道经典小题给定初始高度和反弹恢复系数用代码一步步复现小球的位置变化直到它最终停在水面上。如果你刚接触数值模拟或者想用一段不长的 Python 代码理解“离散采样误差”“事件检测”“停止条件”这些概念那这篇内容正好适合你。我会从最直白的固定步长写法开始再逐步换成更可靠的事件驱动写法中间穿插我实际调试时踩过的几个坑。最后你会发现代码本身不难难的是搞清楚每一步计算到底在近似什么、误差从哪里来。1. 区分两种“乒乓”本文demo要模拟的是哪一种先把这个歧义说清楚否则后面容易对不上号。乒乓球这个比喻在软件里最常见的是“乒乓缓冲”英文叫 ping-pong buffering典型场景是显卡渲染或者数据采集两块缓冲区轮流充当读写目标A 写的时候 B 读下一帧换过来避免读写冲突。这种“换来换去”的结构确实很形象但它是工程架构上的乒乓不是物理上的乒乓。本文的 demo 走的是另一条线即物理里的反弹振荡。一个球从空中被释放受重力往下掉碰到水面后速度反向但幅值会打折扣。因为它每次反弹都损失一部分能量所以反弹高度会越来越低最后无限逼近水面。这个过程在纸上算其实很简单但要用代码模拟就不那么简单了因为你得决定“怎么把连续运动切成离散的时间步”还要处理“小球是否刚好触及水面”这种边界判断。我把目标限定为用代码复现一维竖向上小球的运动轨迹横轴是时间纵轴是高度。初始高度自定反弹恢复系数通常记作 rho自定重力加速度取 9.81。你要观察到的是高度峰值按几何级数递减间隔越来越短最终趋于静止。这个模型虽然简化到只剩一条直线但它已经包含了数值模拟最核心的三件事——离散化、事件判断、终止条件。把这三点弄明白后面你去看更复杂的物理引擎或者游戏跳跃系统思路都是相通的。2. 把一句话问题翻译成物理模型坐标、参数与运动法则很多教程一上来就写代码但我建议先花两分钟把模型定下来不然后面改起来很痛苦。我用 y 表示小球相对水面的高度水面是 y0正方向向上。小球在空中运动时只受重力加速度是 -g。当它碰到水面那一个瞬间速度方向反转并乘以恢复系数 rho也就是 v_new -rho * v_old。rho 取 0.8 表示每次反弹保留 80% 的速度值或者说保留 64% 的动能因为动能和速度平方成正比。这里有个容易忽略的点现实中的反弹并不是瞬时完成的球入水后还有挤压、形变、排水这些过程。但作为理解性 demo我们可以把碰撞当成一个瞬间事件只要是 y0 就立刻反转速度这已经能很好体现“周期衰减”的核心特征了。于是运动法则可以写成下面这几条空中速度以恒定重力加速度递减v v - g * dt位置随速度更新y y v * dt碰到水面y 0把 y 修正回 0然后 v -rho * v你可能想问为什么不能用解析公式直接算其实可以。从高度 h 落到水面的时间是 sqrt(2h/g)反弹后再上升到最高点的时间是 sqrt(2h/g)因为 h rho^2 * h所以下一次触水时间间隔会等比缩小。理论上你甚至能写出等比数列求和公式。但这样做你就失去了“数值模拟”的意义。写代码模拟的价值在于当你加入空气阻力、变加速度、不规则碰撞面之后解析解不再存在你必须依赖逐步推进的计算而这就是数值方法发挥作用的地方。为了防止后面讨论混淆我把代码里用的参数固定下来参数含义取值g重力加速度9.81 m/s^2rho恢复系数0.8y0初始高度1.0 mv0初始速度0.0 m/sdt时间步长0.01 s3. 第一版demo固定时间步长模拟先让它跑起来固定时间步长法是最容易上手的写法。想象你用一个摄像机每隔 dt 秒拍一张照片所有物理量只在快门瞬间被更新。代码里就是一个 while 循环每轮做一次速度更新和一次位置更新最后判断有没有越过水面边界。下面这版是完整的可运行代码import matplotlib.pyplot as plt g 9.81 rho 0.8 dt 0.01 y 1.0 v 0.0 t 0.0 times [t] positions [y] while t 10.0: v v - g * dt y y v * dt t t dt if y 0.0: y 0.0 v -rho * v times.append(t) positions.append(y) plt.plot(times, positions) plt.xlabel(time (s)) plt.ylabel(height (m)) plt.title(bounce motion with fixed dt 0.01) plt.show()跑完这段代码你会看到一条快速衰减的振荡曲线。前几次反弹的高度差非常明显后面则越来越密最终贴近零轴。这不是物理上“停止”而是曲线在图上太密集了看不出变化。控制台里如果把最后几次 v 和 y 打出来你会发现它们还在小数点后好几位继续挣扎。选择 dt0.01 对这类问题已经算比较精细。你可以试着把 dt 调成 0.001曲线变化不大再试着调成 0.1曲线会开始出现肉眼可见的失真。这正是数值模拟里最基础也最核心的训练理解步长对结果的影响。固定步长的好处是逻辑简单所有时间点均匀分布你甚至不用记录每次反弹的具体时刻直接按顺序画图就行。坏处则是碰撞时刻很少能正好落在采样点上于是误差会被反复累积。4. 固定步长跑不动的三处坑穿透、卡死和能量漂移第一版代码能跑但不代表它没问题。我把 dt 从 0.01 改成 0.1 之后马上看到三种诡异的状况这在理解“乒乓原理”时特别容易被忽略值得单独写一节。穿透。小球从最高点下落如果一步跨得太大这一帧 y 还是正的下一帧 y 已经变成 -2.0。你的判断条件 y 0 确实能拦住它但要注意它是在“很深的地方”才被拦住的实际水面位置已经被跳过。这会直接导致反弹高度偏低因为整段负位移对应的势能没被计入。更严重的情况是如果这一步跨出的距离长到让小球的负速度非常大反弹后速度仍然为负下一帧它又立刻撞水面循环看起来就像小球“粘”在水面下方来回颤抖。卡死。恢复系数不是 0 的时候小球反弹速度会越来越小。当 |v| 已经很小而 dt 又不够小时每次更新的位移量 v*dt 可能小于浮点数能分辨的最小单位此时 y 的变化几乎停滞。配合穿透修正“强制把 y 设回 0”小球可能在某几帧内被反复触发反弹逻辑进度却不往前走。表现为曲线尾部出现一段水平抖动而不是平滑收敛。能量漂移。固定步长下碰撞检测永远发生在 y 的采样值小于等于 0 之后这意味着你总是比真实碰撞晚一点才反转速度。这个“晚一点”会让系统多走一段向下的位移等效于多损失了一部分能量。每次反弹都多损失一点整体衰减就比理论上更快。如果你画一条理论峰值衰减线做对照固定步长的那条会明显偏低。要解决这三个问题第一反应是把 dt 调小这当然有效但不能根治因为不管 dt 多小碰撞时刻依旧几乎不可能正好落在采样点上。更稳妥的思路有两个方向一是做“插值补偿”在检测到穿透时用线性插值估算出真实穿透前的位置和速度再在那个近似接触点上处理反弹二是彻底换思路改成事件驱动。5. 来一次事件驱动模拟精确算准每一次触水时刻所谓事件驱动就是不按固定时间间隔走而是先求出“下一次碰撞发生在什么时候”直接跳到那个时刻然后再从碰撞后的状态继续计算。对于纯重力模型这个时间可以用解析式精确算出根本不需要一点一点逼近。假设当前小球在高度 h 处初速度为零或者只有竖直方向速度那么它落到水面所需的时间是 t sqrt(2h/g)触水瞬间的速度是 v -sqrt(2gh)。反弹后的速度是 rho * sqrt(2gh)它能重新上升的最大高度是 h (rho^2 * 2gh) / (2g) rho^2 * h。这样下一次“从最高点落回水面”的时间和这一次完全一样周期规律就出来了。按这个思路代码可以写得很干净import math def event_driven_bounce(y0, rho, g9.81, max_bounces50): times [0.0] heights [y0] h y0 for i in range(max_bounces): if h 1e-6: break t_air math.sqrt(2.0 * h / g) v_impact math.sqrt(2.0 * g * h) t_next times[-1] t_air times.append(t_next) heights.append(0.0) v_rebound rho * v_impact h v_rebound * v_rebound / (2.0 * g) return times, heights times, heights event_driven_bounce(1.0, 0.8) for i in range(min(len(times), 10)): print(fbounce {i}: t {times[i]:.6f}, h {heights[i]:.6f})这个版本的输出会有意思得多。你看到的触水时刻几乎就是教科书等比数列第一次在 0.4515 秒左右第二次在 0.8598 秒左右之后每次间隔都在缩短。如果你把峰值高度提取出来会发现峰值 $h_n$ 正好是 $h_0 * rho^{2n}$。这种精确数字能给 demo 提供很好的验证基准你可以拿它和固定步长法的输出对比立刻就能看出离散化误差有多大。事件驱动版也有局限。它虽然算得准但要求每次碰撞之间的动力学过程有解析解。一旦加入空气阻力或者碰到加速度不是常数的情况你没法再写出 t sqrt(2h/g) 这种式子只能回到固定步长或者混合方法。所以这两套思路并不是谁取代谁而是互补。6. 加入阻尼与停止条件让小球在有限时间内停在水面上纯重力加固定恢复系数的模型在理论上小球永远不会完全静止因为每次反弹高度只是衰减不会变成零只是越来越趋近于零。真实世界当然不是这样因为空气阻力和水的粘性会做功球最终会停在平衡位置。要让 demo 更接近实际可以在运动方程里加一项与速度成正比的阻尼项。阻尼的物理意义很好理解速度越快阻碍越强。代码里加一个阻尼系数 b每次速度更新时写成 v v - (g b * v) * dt。注意这里 b * v 是阻力产生的加速度方向和速度相反所以 b 0 时阻力是让球减速的不管球正在上升还是下落它都在消耗机械能。选 b0.3 时你会发现曲线的衰减速度明显变快后段不再有漫长的“微小振荡”。加上阻尼之后反弹的峰值高度不再是严格的等比数列而是更快地衰减。停止条件也得调整。你不能用 h 1e-6 作为唯一终止条件因为阻尼很强时球可能在高度还比较高的情况下已经几乎没有速度这时候它不会再发生反弹而是直接静止在水面上。所以更合理的停止判断是速度绝对值小于某个阈值且位置离水面足够近比如if abs(v) 0.001 and abs(y) 0.001: break这个阈值怎么选取决于你要模拟的精度。如果你只关心曲线大致形状0.01 就够了如果你要做数值实验对比理论值建议设到 1e-5 以下。我自己的习惯是先把阈值设得宽松跑一遍看大致行为再逐步收窄直到结果不再显著变化这说明阈值已经不影响结论了。7. 从demo到可交差绘图、指标和参数敏感性写完模拟还只是第一步要让这个 demo 真正“能说明问题”最好再加上可视化和几项关键指标。可视化用一种图就够了横轴时间、纵轴高度曲线自然呈现一座座越来越矮的山峰。如果你想更酷一点可以把碰撞点用散点标出来这样事件驱动版里每次触水的时刻就一目了然。下面是给固定步长版加碰撞点标记的代码片段hits_t [] hits_y [] # 循环内检测到 y 穿越水面时记录 if y 0: hits_t.append(t) hits_y.append(0.0)指标方面建议至少计算两个值。第一个是第一次反弹峰值它能直接验证 rho 的准确性第一次峰值应接近 rho^2 * y0如果偏差超过 5%说明步长太大或阻尼加得过猛。第二个是相邻碰撞的时间间隔在无阻尼事件驱动模型里间隔比例应该是 rho如果这个比例不对说明碰撞判断不精确。参数敏感性实验是让这个 demo 变得“有深度”的最快方法。你可以固定 g 和 y0把 rho 从 0.5 扫到 0.95画出若干条曲线叠在一起。这时你会发现 rho 越大曲线尾部拖得越长振荡次数越多这和直觉一致恢复系数越接近 1能量损耗越少球越难停下来。再把 dt 分别设为 0.1、0.01、0.001 跑同一组参数曲线之间的差异就是数值误差的可视化证据。这些图做出来后这个 demo 就从一个“能跑的代码”变成了一个“能讲清楚的演示”。面试、汇报或者写作业的时候光贴代码是没有说服力的但三张图并排一摆误差、参数影响、收敛行为全都清楚了。8. 基于反弹模型还能玩出哪些花样以及我的收尾建议这套乒乓反弹的核心逻辑其实可以向外延伸很多。比如加上双球两个小球同时下落一个恢复系数 0.6一个 0.9初始高度相同观察它们什么时候反弹会交错。比如改成侧向平面初始速度带水平分量水面换成地面小球变成弹跳球这已经接近游戏里角色跳跃的物理原型。再比如把水面换成移动的活塞面球撞上去后的反弹速度还要叠加活塞的运动速度这就变成了一个更接近工程实际的碰撞问题。代码层面如果你用的是 Python还可以尝试把固定步长和事件驱动两种方式封装成两个函数输入输出统一方便直接对比。这种“先做对再做好”的思路比一开始就追求复杂模型更有价值。我自己的实践里最常见的错误就是过早加入各种花哨的外力导致失败时根本定位不到是碰撞逻辑的问题还是数值发散的