尧图网站建设 
这几年我一直在用Python折腾物理模拟从最初做几个小球弹跳的网页小玩具到后来给机器人控制算法搭动力学验证环境再到给虚拟现实培训项目做碰撞反馈原型绕了一大圈最趁手的工具还是Python加NumPy。很多同行一听到“用Python写物理引擎”第一反应就是“慢”“玩具”“不专业”。但如果你认真算过一笔账——开发周期、迭代效率、和机器学习生态的对接便利性——就会同意我的判断NumPy把Python的数值计算性能拉高了两到三个数量级之后很多原本以为必须上C的活Python其实完全扛得下来。这篇文章会把我从零搭一套物理模拟引擎的全过程摊开讲从环境配置、核心数据结构设计到粒子系统、刚体动力学、碰撞检测与约束求解再到性能优化和一个能跑的完整沙盒。游戏开发、机器人控制、虚拟现实这些方向的朋友尤其是想让手头仿真代码跑得更顺的同学应该都能从这里抄到一些作业。1. 为什么选PythonNumPy做物理引擎先解决“慢”的成见1.1 物理引擎到底在算什么三大核心模块先说清楚物理引擎的本质。无论多炫酷的引擎每一帧做的事情翻来覆去就是三件动力学积分、碰撞检测、约束求解。动力学积分是让物体动起来本质就是对着牛顿第二定律做数值积分已知质量、受力算加速度然后更新速度和位置。碰撞检测是判断“谁撞了谁”从粗略的空间剔除Broad Phase到精确的形状相交测试Narrow Phase每一步都在做距离、投影、相交这类几何运算。约束求解是让结果变得真实两个球不能互相嵌入物体不能穿墙铰链和关节得保持连接关系——这些约束通常会转成一堆线性方程组或者迭代修正问题。这三个模块放到Python和NumPy的视角下看全部都是数组运算。位置是一批(N,2)的数组速度是(N,2)的数组距离检测是一个广播减法加一次按行求和约束修正是一次向量化的乘法加法。这正是NumPy的主场也是整个技术方案成立的根本原因。1.2 什么项目真正适合用Python写这里我不想吹“Python万能”因为它确实不是。我按自己的经验给一个判断标准当单场景物体数量在几万以内、物理步长能跑到120Hz、且需要频繁调整算法或对接AI训练流程时PythonNumPy的性价比碾压C当目标是落地到AAA游戏或者超大开放世界时才需要认真考虑原生引擎。我画过一张很粗的判断表直接抄在这里场景物体规模适合方案原因教学演示、算法验证小于5000PythonNumPy改造成本极低随时改参数跑结果机器人控制策略训练小于1000PythonNumPy直接嵌入强化学习环境天然无缝VR培训原型、体感交互小于3000PythonNumPy迭代快方便调手感大型游戏落地、商业引擎数万以上C/Rust内存和指令级优化要求极高这背后的逻辑是仿真开发里真正烧时间的往往不是运行而是“调参、看现象、发现问题、再改代码”这个循环。Python把循环的周期从几个小时压缩到几分钟这个优势比单帧运行快几十毫秒重要得多。1.3 NumPy为什么比普通列表快这么多热词里有“numpy和list比快在哪”这个问题我正好在一次分享里实测过。用一百万元素做逐元素加法普通Python列表用列表推导式跑一遍NumPy数组用加法运算符跑一遍计时结果通常是列表0.08秒左右NumPy 0.001秒出头差距五十倍上下。这个差距的根源有三层普通列表里存的是指向Python对象的指针每个对象有类型信息、引用计数循环时Python解释器要一行一行地把字节码翻译成操作整条链路非常重。NumPy数组则把数据连续地铺在内存里底层是C写的循环编译器还能自动向量化用上CPU的SIMD指令一次性批处理多个数据。第三层是缓存友好性连续内存访问对CPU缓存极友好而链表式的对象散布会把缓存命中率打到很低。做物理引擎时我们恰恰要把所有对象的状态连续存放然后对整批数据做同样的运算——这几乎是为NumPy量身定做的模式。2. 环境配置与数据结构开工前的关键准备2.1 把开发环境一次配好含常见安装报错很多卡在第一步的人问题几乎都出在“多个Python环境混着用”。我的建议是不要直接把numpy装到系统Python里用一个独立的虚拟环境隔离项目依赖。从零开始的操作顺序是下载安装Python装的时候注意勾选“Add Python to PATH”然后打开终端创建并激活虚拟环境python -m venv .venv # Windows .venv\Scripts\activate # macOS / Linux source .venv/bin/activate激活后安装依赖pip install --upgrade pip pip install numpy几个高频报错我顺手列出来都是热词里出现率极高的ModuleNotFoundError: No module named numpy大概率是解释器没切对。VSCode按下快捷键CtrlShiftP输入“Python: Select Interpreter”选中刚才创建的虚拟环境。PyCharm里在File Settings Project Python Interpreter里把解释器换成.venv路径。安装numpy时卡在“Installing backend dependencies”这是老版本pip在尝试走源码编译构建后端。解决办法是先把pip升级到最新再装或者指定只装预编译二进制包pip install --only-binary :all: numpy。绝大多数情况装完之后立刻就好。numpy版本不匹配典型场景是新版Python配了旧版Numpy。比如Python 3.13配上numpy 1.24就会报警告甚至直接报错。无脑做法是pip install --upgrade numpy或者到PyPI看当前最新版本号直接指定安装。还有个小建议如果打算长期做仿真开发把环境写好到一个requirements.txt里numpy1.26 numba0.59 matplotlib下次换机器执行pip install -r requirements.txt就完事了。2.2 用数组而不是对象来管理物理世界这是整个引擎设计里最值得讲清楚的一步。新手很容易把每个小球或者刚体写成一个class每个物体持有自己的位置、速度然后在一个对象列表里循环。这个写法直观但性能极差而且越改越乱。我采用的方式是平行的数组结构——每个物理量是一个独立的NumPy数组数组的第i行就代表第i个物体import numpy as np N 1000 pos np.zeros((N, 2), dtypenp.float64) # 位置 vel np.zeros((N, 2), dtypenp.float64) # 速度 force np.zeros((N, 2), dtypenp.float64) # 合力 inv_mass np.ones(N) # 逆质量 radius np.full(N, 0.05) # 半径这里有个很重要的细节用逆质量inv_mass而不是质量mass。动力学公式里加速度 力除以质量写成逆质量就是力乘以inv_mass。为什么不用质量因为逆质量等于0时表示“质量无穷大的静态物体”比如地面、墙壁它永远不移动。如果存的是质量就得在代码里到处判断“质量是否为零”而用逆质量后乘法直接天然处理零乘任何数都是零干净利落。用这种结构做整批更新代码会非常优雅。施加重力并更新速度只需要五行acc force * inv_mass[:, None] vel acc * dt pos vel * dt force.fill(0.0)inv_mass[:, None]是把(N,)形状扩展成(N,1)让它可以和(N,2)的force逐元素相乘。这类形状操作是NumPy编程的必修手感。2.3 引擎里常用的线性代数操作物理引擎逃不开矩阵运算。这里列几个我在代码里真正用过的不是抄文档矩阵乘法用符号。刚体的旋转矩阵更新、刚体局部点到世界坐标的变换全部用R local_point完成。二维刚体旋转矩阵是2x2和(N,2)批量点做乘积时pos_world pos_local R.T。np.linalg.inv解约束方程组、np.linalg.det判断几何退化。比如在分析场景中计算平面三角形的面积行列式就是面积的二倍。np.einsum做批量内积这个我强烈推荐。np.einsum(ij,ij-i, d, d)可以一次算出每一行向量的平方长度等效于np.sum(d*d, axis1)但可读性更强。热词里的“numpy三维数组相乘”也适用当你要并行模拟K个独立的物理世界时状态可以组织成(K, N, D)三维数组用np.matmul或者做批量矩阵乘法一次更新所有世界。这在强化学习里批量采样训练数据时特别常用。3. 第一个物理模块粒子系统与稳定积分3.1 十几行代码让500个小球下落粒子系统是物理引擎的最小单元拿来验证NumPy的写法和积分的稳定性再合适不过。下面是我跑通的第一版代码就十几行核心逻辑import numpy as np N 500 dt 1.0 / 120.0 g np.array([0.0, -9.8]) pos np.random.uniform([0.0, 0.0], [4.0, 3.0], (N, 2)) vel np.zeros((N, 2)) radius np.full(N, 0.03) for frame in range(240): # 重力加速度 vel g * dt # 位置更新 pos vel * dt # 地面碰撞反弹 floor_y 0.0 below pos[:, 1] floor_y radius pos[below, 1] 2 * (floor_y radius[below]) - pos[below, 1] vel[below, 1] -0.8 * vel[below, 1]注意看第14行到第16行的地面处理这里用了一个镜像反弹技巧求越界位置相对于边界的镜像位置而不是简单地把y坐标塞回边界值。这样能保证速度和位置的关系不产生新的能量突变。这段代码已经是完整的“粒子重力边界碰撞”物理系统用matplotlib.animation加个画图就能直接看到500个球在下落反弹。3.2 三种积分器怎么选数值积分的选择直接决定整个引擎会不会“爆”。我用过三种把各自特性说清楚显式欧拉是先更新位置再用新位置算加速度或者先算位置再更新速度两种变体在有弹簧这类力时会持续向系统注入能量模拟几十帧后整个场景就开始发疯。半隐式欧拉也叫symplectic Euler是先更新速度、再用新速度更新位置vel acc * dt pos vel * dt这个顺序的小改动带来质变因为对于自由落体和弹簧系统半隐式欧拉能保持系统能量上下波动而不是单调发散。这是游戏物理默认选择也是我所有后续代码的基础。Verlet积分用前一帧位置和当前帧位置做差分new_pos 2 * pos - old_pos acc * dt * dt它没有一个显式的速度量但天然适合约束很多的场景——布料、绳索、粒子堆叠。因为约束求解时只调整位置、从位置变化反推速度误差累积特性比直接操作速度好很多。选型建议初学阶段无脑用半隐式欧拉加固定时间步长。固定步长很关键可变步长会让物理系统在不同分辨率下表现不一致这是很多诡异抖动的根源。3.3 墙壁反弹与恢复系数边界处理除了镜像反弹还要引入恢复系数e它描述碰撞后法向速度保留多少。e1是完全弹性碰撞能量不损失e0是碰撞后法向速度清零物体贴住地面不动。实际场景一般取0.5到0.9之间代码里就是一个系数vel[below, 1] -0.8 * vel[below, 1]有个很容易被忽略的问题当球的落速极小时反弹速度也极小物体停不下来会一直在地面附近蠕动视觉上很烦。解决方法是给一个休眠阈值速度低于某个值直接清零threshold 0.05 low_speed np.abs(vel[below, 1]) threshold vel[low_speed, 1] 0.0这个细节在游戏手感里特别重要我后来做交互原型时几乎每个项目都用到。4. 刚体动力学与碰撞约束从“点”到“块”4.1 刚体的状态多出来的旋转自由度粒子只有两个自由度x和y刚体多一个旋转自由度θ。二维刚体的状态量变成五个位置(x, y)、速度(vx, vy)、角度θ、角速度ω。同样用数组组织pos np.zeros((N, 2)) vel np.zeros((N, 2)) angle np.zeros(N) # 物体朝向 omega np.zeros(N) # 角速度 I np.zeros(N) # 转动惯量转动惯量I是旋转动力学里的“质量”。实心圆盘对质心轴的转动惯量是0.5 * m * r^2它决定了同样力矩下角加速度有多大。刚体局部坐标系里的一个点r_local世界坐标计算公式是R np.array([[np.cos(angle[i]), -np.sin(angle[i])], [np.sin(angle[i]), np.cos(angle[i])]]) world_point pos[i] R r_local如果有一个力F作用在物体上但作用点不在质心它除了产生线加速度还会产生力矩τ。二维里的力矩是一个标量τ rx * Fy - ry * Fx其中(rx, ry)是从质心指向受力点的向量。这个公式在实现“推力器”“风帆”这类玩法时非常有用。4.2 碰撞检测Broad Phase 与 Narrow Phase碰撞检测是物理引擎最容易拖慢性能的部分。朴素做法是遍历所有物体两两比较复杂度O(N^2)500个球每帧要做12.5万次距离判断5000个球就是1250万次完全不可持续。工业级做法分成两阶段。Broad Phase快速排除明显不相交的物体对把候选对数量降到最小。最简单实用的是空间网格格子哈希把空间划分成边长略大于最大物体尺寸的均匀格子每个物体登记到自己所在的格子然后只检查同格子和相邻九格内的物体。因为物体密度通常不高平均复杂度能降到接近O(N)。Narrow Phase再精确判断候选对是否真的相交。球和球用距离平方和比较凸多边形用分离轴定理SAT——依次检查每条边的法线方向如果存在任何一个方向上两个图形的投影没有重叠就说明图形分离。这个定理是不相交判定逻辑很干净所有凸多边形碰撞的基础都是它。我在实际项目里的建议是先用朴素O(N^2)把逻辑跑通再用网格优化。过早优化会让代码复杂度翻倍掩盖了物理逻辑本身的问题。4.3 约束求解冲量、位置修正和一招稳定术碰撞发生后要做两件事阻止互相穿透、改变运动方向。经典做法是冲量法。两个球碰撞时沿法线方向的相对接近速度为vn碰撞后要变成方向相反、大小为e倍。解这一个方程得到法向冲量的公式j -(1 e) * vn / (inv_mass_i inv_mass_j)然后把冲量施加到两个物体上vel[i] j * inv_mass_i * n vel[j] - j * inv_mass_j * n这个公式在等质量弹性碰撞下正好退化成“交换法向速度”物理直觉完全对得上。但光改速度还不够离散时间步会让两个球在检测出碰撞时已经嵌入了一小段距离。只看速度不修正位置时间一长物体就会“沉”到地面里。所以我加了位置修正把重叠量按质量比例分摊给两个物体各推回一半质量大的动得少质量小的动得多。交给读者一个小技巧约束解算不要只跑一遍。很多引擎同一帧内会迭代2到4次约束修正每次把残留误差减半。迭代3次的效果和迭代20次接近而性能开销小得多。这也是为什么很多“堆叠箱子”演示里底部箱子看起来稳如泰山——其实是迭代约束的功劳。再进一步就是位置动力学PBD先根据力预测一个新位置然后把所有约束当成投影操作直接在位置上修正最后用修正后的位置反推速度。这套思路写布料、绳索、果冻这类软体模拟几乎是无脑好用的。4.4 提高模拟稳定性的几个工程细节经验的坑踩多了之后我习惯性地在刚体系统里加这些保护最大速度限制。每帧更新完速度后把所有速度幅值钳制在一个上限内比如20 m/s。这能防止离散碰撞检测漏检高速物体也能阻止数值溢出。阻尼项。每帧乘一个略小于1的系数比如vel * 0.999。它模拟空气阻力同时把数值误差缓慢压掉。休眠机制。速度低于阈值且连续数帧不变就把物体标记为休眠跳过它的积分和碰撞计算。大量堆叠刚体稳定后休眠能省下90%以上的计算量。物体的摩擦力近似。在碰撞速度修正里把切向速度乘以一个摩擦系数再保留一部分能产生粗糙表面的效果不需要引入复杂的摩擦锥模型。这些在引擎代码里往往只有几行但对结果的影响是决定性的。5. 性能优化向量化、预分配与Numba5.1 消灭for循环的两条路线Python物理引擎的性能瓶颈九成在for循环。消灭它们有两条路线。第一条是批量向量化——把外层循环去掉一次对整批数据做运算。比如施加外力代码从for i in range(N): force[i] ...变成force ...整个数组一次性算。或者把多个物体的运动统一写成矩阵运算一次执行。第二条是**“半向量化”**外层保留一个轻量循环循环内全部用向量操作。以球-球碰撞检测为例对于每个球i一次计算它和其他所有球的距离矩阵d pos - pos[i] # 广播所有球相对球i的偏移 dist_sq np.einsum(ij,ij-i, d, d) # 每个球到球i的距离平方 min_d radius[i] radius collide (dist_sq 0) (dist_sq min_d * min_d) idx np.flatnonzero(collide) # 找到真正碰撞的球索引这样外层循环N次内层全是C级别速度的向量操作比纯Python的两重循环快一到两个数量级而且可读性依然很好。5.2 减少内存分配的三个习惯NumPy数组加法a b会生成一个全新的临时数组物理引擎每帧要做大量这样的操作临时数组多了内存分配和释放本身就会拖慢速度。三个习惯能明显改善第一能原地操作就原地操作。vel g * dt这种是原地操作不产生新数组但如果写vel vel g * dt就产生新对象。第二带out参数计算。np.add(a, b, outc)把结果写入预分配的c省掉中间临时量。第三避免在循环里拼接数组。反复调用np.concatenate会让内存不断搬迁正确做法是先分配最大容量再对有效部分切片赋值。这三个习惯都改完之后我之前一个波形仿真项目性能提升了30%左右代码逻辑完全没变。5.3 Numba在不改物理代码的前提下再快一个数量级如果需要更大提升我推荐Numba——一个直接把Python函数编译成机器码的JIT库。它对NumPy的支持相当好很多Python循环代码只要加个装饰器就能提速几十倍。把积分和碰撞检测写成一个函数加上njitfrom numba import njit njit(cacheTrue) def physics_step(pos, vel, radius, inv_mass, dt): N pos.shape[0] vel g * dt # 球与球碰撞的纯Python循环 for i in range(N): for j in range(i 1, N): dx pos[i, 0] - pos[j, 0] dy pos[i, 1] - pos[j, 1] dist_sq dx * dx dy * dy min_d radius[i] radius[j] if dist_sq min_d * min_d and dist_sq 0: # 速度冲量计算 ... pos vel * dt这里有意思的地方在于Numba允许你写回“朴素的二重循环”因为编译后根本不再经过Python解释器。这比写一堆精妙的NumPy向量化技巧还快并且代码更直白。我自己处理几千个刚体碰撞时纯Python版跑2000帧要十几秒改成Numba后两秒以内稳跑。5.4 用一张表看三种实现方式的差距拿“600个小球堆积在箱子里模拟1000帧”做了一个很粗糙的基准三类写法的差距大约是实现方式骨架难度1000帧耗时我的旧笔记本可读性纯Python二重循环低22秒左右最好NumPy半向量化中3.5秒左右良好Numba JIT编译中0.8秒左右良好优化后的C高0.3秒左右较差Numba方案在代码量接近Python的前提下性能已经摸到了C的脚后跟。对于大多数教学、机器人训练和VR原型场景这个速度足够了。这套优化路径是我个人最喜欢的方向先追求逻辑清晰再考虑向量化最后上JIT每一步都有明确的收益。6. 实战500个球的可交互物理沙盒6.1 场景设计与参数设定理论讲完上一段真正能跑的代码。我设计了一个2D物理沙盒500个半径不等的球落入一个封闭矩形空间受重力影响彼此碰撞、和墙壁碰撞并反弹。关键参数取这些值时间步长1/120秒这样每帧画面有约8个物理子步视觉平滑且较稳定恢复系数0.85模拟出略带弹性的碰撞世界尺寸4米乘3米配合v0.05左右的球半径一个画面里能塞下足够的物体量。为什么是120Hz而不是60Hz因为物理系统对时间步长的敏感度和画面刷新率不是一回事。物理子步越小离散检测穿透的风险越低约束的收敛性越好。游戏引擎里把物理频率定在120Hz或240Hz很常见哪怕画面只刷新60次。6.2 完整代码实现我把它整理成一个完整可运行的脚本import numpy as np # 参数 N 500 DT 1.0 / 120.0 W, H 4.0, 3.0 G np.array([0.0, -9.8]) E 0.85 DAMP 0.999 # 状态 pos np.random.uniform([0.1, 0.1], [W - 0.1, H - 0.1], (N, 2)).astype(np.float64) vel np.zeros((N, 2)) radius np.random.uniform(0.03, 0.06, N) inv_mass np.ones(N) def step(): global pos, vel # 重力 空气阻尼 vel G * DT vel * DAMP # 球与球碰撞 for i in range(N): d pos - pos[i] # (N, 2) dist_sq np.einsum(ij,ij-i, d, d) min_d radius[i] radius collide (dist_sq 0) (dist_sq min_d * min_d) idx np.flatnonzero(collide) if idx.size 0: continue dist np.sqrt(dist_sq[idx]) n_ij d[idx] / dist[:, None] # 从i指向j的单位向量 rel_v vel[i] - vel[idx] vn np.sum(rel_v * n_ij, axis1) approaching vn 0 # 正在接近才处理 idx_app idx[approaching] if idx_app.size 0: continue n_app n_ij[approaching] vn_app vn[approaching] # 冲量 denom inv_mass[i] inv_mass[idx_app] j -(1 E) * vn_app / denom vel[i] np.sum((j * inv_mass[i])[:, None] * n_app, axis0) vel[idx_app] - (j * inv_mass[idx_app])[:, None] * n_app # 位置修正按质量比分摊重叠 overlap (radius[i] radius[idx_app]) - dist[approaching] w_i inv_mass[i] w_j inv_mass[idx_app] total_inv w_i w_j if total_inv 0: pos[i] np.sum((-overlap * w_i / total_inv)[:, None] * n_app, axis0) pos[idx_app] (overlap * w_j / total_inv)[:, None] * n_app # 推进位置 pos vel * DT # 墙壁约束与反弹 left pos[:, 0] radius right pos[:, 0] W - radius bottom pos[:, 1] radius top pos[:, 1] H - radius pos[left, 0] 2 * radius[left] - pos[left, 0] vel[left, 0] -E * vel[left, 0] pos[right, 0] 2 * (W - radius[right]) - pos[right, 0] vel[right, 0] -E * vel[right, 0] pos[bottom, 1] 2 * radius[bottom] - pos[bottom, 1] vel[bottom, 1] -E * vel[bottom, 1] pos[top, 1] 2 * (H - radius[top]) - pos[top, 1] vel[top, 1] -E * vel[top, 1] # 低速休眠阈值 vel[np.abs(vel) 0.001] 0.0 # 跑1200帧 for frame in range(1200): step()在Notebook里配合matplotlib.animation每50帧画一张散点图就能看到几百个球翻滚落稳的全过程。6.3 关键代码逐段讲解碰撞处理是整段代码的难点我拆开讲几个重点单位向量与接近判断。d pos - pos[i]计算的是“从球i指向其他球”的向量除以距离后就得到每个碰撞对象指向当前球i的单位向量n_ij。rel_v vel[i] - vel[idx]是球i相对碰撞对象的速度把它投影到n_ij上得到法向相对速度vn。只有vn大于0才说明两个球正在接近小于0说明正在分离无需处理——这个判断省掉了大量无用的速度修正。冲量施加。j -(1 E) * vn_app / denom计算出的是一个标量冲量沿着法线方向。球i收到的速度变化是j乘以自身的逆质量再乘方向向量碰撞对象则相反减去同样的量。因为这里是批量计算用了[:, None]扩展维度配合广播。位置修正的权重分配。重叠量必须消除但让谁动得多、谁动得少按逆质量比例分配。静态物体inv_mass为0时total_inv可能为0需要判断否则会除零。这个判断虽然小但少了它一旦混入静态物体引擎会瞬间崩掉。墙壁镜像反弹。墙壁部分没有用循环直接对布尔掩码索引的数组做赋值这正是NumPy风格的体现。镜像公式在边界之外也能正确处理不会卡死。6.4 向游戏、机器人、VR场景扩展沙盒能跑之后扩展方向已经很明朗。游戏方向把多边形的碰撞检测补上SAT加入关节约束旋转铰链、滑橇就可以做平台跳跃的谜题原型或者模拟弹弓、荡绳这些小玩法。Unity和Godot的物理引擎当然比这个完善但自己写一遍才能理解引擎参数背后的含义调起手感来也更准。机器人控制方向把动力学部分替换成机械臂的刚体链加上关节力矩就能和强化学习环境对接。我用Gymnasium自定义环境时核心的step函数里就是一套和上面结构完全相同的NumPy物理更新。批量状态处理用3D数组还能并行开多个环境的训练采样效率很高。VR培训方向沙盒里的碰撞反馈逻辑可以直接转成手势交互原型。比如模拟工具碰触虚拟物体时的震动强度用碰撞冲量大小作为反馈源的参数用户体验会非常自然。7. 常见问题排查与避坑实录7.1 环境类问题的速查表我把开箱时碰到的高频问题整理成表方便以后直接翻问题可能原因解决办法No module named numpy当前解释器不是虚拟环境VSCode里切换解释器或重新激活.venv安装卡在backend dependenciespip版本太旧尝试源码构建升级pip或--only-binary :all:强制用预编译包numpy版本不匹配报错Python版本太新、numpy太旧pip install --upgrade numpy代码在A电脑正常B电脑报错依赖版本不一致用requirements.txt锁定版本号运行中发现结果不确定用了随机数却没固定种子初始化时np.random.seed(42)7.2 物理模拟数值不稳定的几个典型症状症状一堆叠物体持续抖动、缓慢漂移。八成是约束迭代次数不够或时间步长太大。先把物理步长降到120Hz再把碰撞位置修正迭代次数提到4你会发现整个场景像被按住了。症状二快速飞行的球直接穿透墙壁。这是离散碰撞检测的经典缺陷——物体在一个时间步内从墙壁一侧运动到另一侧检测完全没碰到。解法有收紧时间步、给物体加连续碰撞检测、或者限制最大速度。我常用后两者组合。症状三整体能量越变越大物体会自己飘起来。这是积分器不稳定加碰撞处理给系统注入能量。检查是不是用了纯显式欧拉改成半隐式检查恢复系数是不是设成了大于1最后加上阻尼项兜底。症状四刚体旋转起来越来越快最后失控。通常是转动惯量算错或者力矩计算时把作用点向量用反了。二维里力矩τrxFy-ryFx方向正负搞反会让角速度往错的方向累加。7.3 我踩过的最贵的几个坑第一个坑是用逆质量表格处理约束时忘掉静态物体。第一次让箱子堆到地面后一个除零直接把整个进程送走。后来养成了习惯任何涉及inv_mass的除法先判断分母。第二个坑是为了追求全向量化把碰撞对全部构造成pair矩阵。1000个球会产生50万行pair中间矩阵瞬间吃掉几个G内存直接被系统警告。这是个很典型的“向量化思维过猛”案例真实场景应该用broad phase配合半向量化而不是把NPair全铺开。第三个坑是把物理步长和渲染帧率绑定。早期代码直接每帧调一次step结果在低帧率机器上物理明显变慢在高帧率机器上又太快。后来统一改成固定物理步长加累计器把渲染和物理完全解耦所有屏幕上的表现才一致起来。第四个坑比较深刻把所有物体集中到少量大数组之后调试几何问题变得极其痛苦。某个球的“位置”在数组里只是其中一行出了诡异行为没法直观看到。后来我自己习惯是核心计算用数组但排查问题时写一个def debug_obj(i): print(pos[i], vel[i], radius[i])之类的辅助函数或者临时只保留两三个物体跑最小复现场景。大多数情况下把问题缩小到最少重现条件原因一眼就能看出来。我个人在这些项目里的体会是物理引擎最有意思的不是把公式写对而是在“写对公式之后”的那个晚上——你盯着屏幕上几百个球自己推来推去发现它们有了沉甸甸的质感。PythonNumPy这个组合恰好能让你在最短时间内走到那个晚上。先把这套基础跑通后面不管是加关节、加流体还是加软体都是顺着同一套数组思维延伸出去的事。关于版本兼容性的问题建议换新机器时先跑一遍沙盒脚本确认环境一致再继续开发能省下一整天的排查时间。 用Python和NumPy从零搭建物理模拟引擎：碰撞检测、约束求解与性能优化
这几年我一直在用Python折腾物理模拟从最初做几个小球弹跳的网页小玩具到后来给机器人控制算法搭动力学验证环境再到给虚拟现实培训项目做碰撞反馈原型绕了一大圈最趁手的工具还是Python加NumPy。很多同行一听到“用Python写物理引擎”第一反应就是“慢”“玩具”“不专业”。但如果你认真算过一笔账——开发周期、迭代效率、和机器学习生态的对接便利性——就会同意我的判断NumPy把Python的数值计算性能拉高了两到三个数量级之后很多原本以为必须上C的活Python其实完全扛得下来。这篇文章会把我从零搭一套物理模拟引擎的全过程摊开讲从环境配置、核心数据结构设计到粒子系统、刚体动力学、碰撞检测与约束求解再到性能优化和一个能跑的完整沙盒。游戏开发、机器人控制、虚拟现实这些方向的朋友尤其是想让手头仿真代码跑得更顺的同学应该都能从这里抄到一些作业。1. 为什么选PythonNumPy做物理引擎先解决“慢”的成见1.1 物理引擎到底在算什么三大核心模块先说清楚物理引擎的本质。无论多炫酷的引擎每一帧做的事情翻来覆去就是三件动力学积分、碰撞检测、约束求解。动力学积分是让物体动起来本质就是对着牛顿第二定律做数值积分已知质量、受力算加速度然后更新速度和位置。碰撞检测是判断“谁撞了谁”从粗略的空间剔除Broad Phase到精确的形状相交测试Narrow Phase每一步都在做距离、投影、相交这类几何运算。约束求解是让结果变得真实两个球不能互相嵌入物体不能穿墙铰链和关节得保持连接关系——这些约束通常会转成一堆线性方程组或者迭代修正问题。这三个模块放到Python和NumPy的视角下看全部都是数组运算。位置是一批(N,2)的数组速度是(N,2)的数组距离检测是一个广播减法加一次按行求和约束修正是一次向量化的乘法加法。这正是NumPy的主场也是整个技术方案成立的根本原因。1.2 什么项目真正适合用Python写这里我不想吹“Python万能”因为它确实不是。我按自己的经验给一个判断标准当单场景物体数量在几万以内、物理步长能跑到120Hz、且需要频繁调整算法或对接AI训练流程时PythonNumPy的性价比碾压C当目标是落地到AAA游戏或者超大开放世界时才需要认真考虑原生引擎。我画过一张很粗的判断表直接抄在这里场景物体规模适合方案原因教学演示、算法验证小于5000PythonNumPy改造成本极低随时改参数跑结果机器人控制策略训练小于1000PythonNumPy直接嵌入强化学习环境天然无缝VR培训原型、体感交互小于3000PythonNumPy迭代快方便调手感大型游戏落地、商业引擎数万以上C/Rust内存和指令级优化要求极高这背后的逻辑是仿真开发里真正烧时间的往往不是运行而是“调参、看现象、发现问题、再改代码”这个循环。Python把循环的周期从几个小时压缩到几分钟这个优势比单帧运行快几十毫秒重要得多。1.3 NumPy为什么比普通列表快这么多热词里有“numpy和list比快在哪”这个问题我正好在一次分享里实测过。用一百万元素做逐元素加法普通Python列表用列表推导式跑一遍NumPy数组用加法运算符跑一遍计时结果通常是列表0.08秒左右NumPy 0.001秒出头差距五十倍上下。这个差距的根源有三层普通列表里存的是指向Python对象的指针每个对象有类型信息、引用计数循环时Python解释器要一行一行地把字节码翻译成操作整条链路非常重。NumPy数组则把数据连续地铺在内存里底层是C写的循环编译器还能自动向量化用上CPU的SIMD指令一次性批处理多个数据。第三层是缓存友好性连续内存访问对CPU缓存极友好而链表式的对象散布会把缓存命中率打到很低。做物理引擎时我们恰恰要把所有对象的状态连续存放然后对整批数据做同样的运算——这几乎是为NumPy量身定做的模式。2. 环境配置与数据结构开工前的关键准备2.1 把开发环境一次配好含常见安装报错很多卡在第一步的人问题几乎都出在“多个Python环境混着用”。我的建议是不要直接把numpy装到系统Python里用一个独立的虚拟环境隔离项目依赖。从零开始的操作顺序是下载安装Python装的时候注意勾选“Add Python to PATH”然后打开终端创建并激活虚拟环境python -m venv .venv # Windows .venv\Scripts\activate # macOS / Linux source .venv/bin/activate激活后安装依赖pip install --upgrade pip pip install numpy几个高频报错我顺手列出来都是热词里出现率极高的ModuleNotFoundError: No module named numpy大概率是解释器没切对。VSCode按下快捷键CtrlShiftP输入“Python: Select Interpreter”选中刚才创建的虚拟环境。PyCharm里在File Settings Project Python Interpreter里把解释器换成.venv路径。安装numpy时卡在“Installing backend dependencies”这是老版本pip在尝试走源码编译构建后端。解决办法是先把pip升级到最新再装或者指定只装预编译二进制包pip install --only-binary :all: numpy。绝大多数情况装完之后立刻就好。numpy版本不匹配典型场景是新版Python配了旧版Numpy。比如Python 3.13配上numpy 1.24就会报警告甚至直接报错。无脑做法是pip install --upgrade numpy或者到PyPI看当前最新版本号直接指定安装。还有个小建议如果打算长期做仿真开发把环境写好到一个requirements.txt里numpy1.26 numba0.59 matplotlib下次换机器执行pip install -r requirements.txt就完事了。2.2 用数组而不是对象来管理物理世界这是整个引擎设计里最值得讲清楚的一步。新手很容易把每个小球或者刚体写成一个class每个物体持有自己的位置、速度然后在一个对象列表里循环。这个写法直观但性能极差而且越改越乱。我采用的方式是平行的数组结构——每个物理量是一个独立的NumPy数组数组的第i行就代表第i个物体import numpy as np N 1000 pos np.zeros((N, 2), dtypenp.float64) # 位置 vel np.zeros((N, 2), dtypenp.float64) # 速度 force np.zeros((N, 2), dtypenp.float64) # 合力 inv_mass np.ones(N) # 逆质量 radius np.full(N, 0.05) # 半径这里有个很重要的细节用逆质量inv_mass而不是质量mass。动力学公式里加速度 力除以质量写成逆质量就是力乘以inv_mass。为什么不用质量因为逆质量等于0时表示“质量无穷大的静态物体”比如地面、墙壁它永远不移动。如果存的是质量就得在代码里到处判断“质量是否为零”而用逆质量后乘法直接天然处理零乘任何数都是零干净利落。用这种结构做整批更新代码会非常优雅。施加重力并更新速度只需要五行acc force * inv_mass[:, None] vel acc * dt pos vel * dt force.fill(0.0)inv_mass[:, None]是把(N,)形状扩展成(N,1)让它可以和(N,2)的force逐元素相乘。这类形状操作是NumPy编程的必修手感。2.3 引擎里常用的线性代数操作物理引擎逃不开矩阵运算。这里列几个我在代码里真正用过的不是抄文档矩阵乘法用符号。刚体的旋转矩阵更新、刚体局部点到世界坐标的变换全部用R local_point完成。二维刚体旋转矩阵是2x2和(N,2)批量点做乘积时pos_world pos_local R.T。np.linalg.inv解约束方程组、np.linalg.det判断几何退化。比如在分析场景中计算平面三角形的面积行列式就是面积的二倍。np.einsum做批量内积这个我强烈推荐。np.einsum(ij,ij-i, d, d)可以一次算出每一行向量的平方长度等效于np.sum(d*d, axis1)但可读性更强。热词里的“numpy三维数组相乘”也适用当你要并行模拟K个独立的物理世界时状态可以组织成(K, N, D)三维数组用np.matmul或者做批量矩阵乘法一次更新所有世界。这在强化学习里批量采样训练数据时特别常用。3. 第一个物理模块粒子系统与稳定积分3.1 十几行代码让500个小球下落粒子系统是物理引擎的最小单元拿来验证NumPy的写法和积分的稳定性再合适不过。下面是我跑通的第一版代码就十几行核心逻辑import numpy as np N 500 dt 1.0 / 120.0 g np.array([0.0, -9.8]) pos np.random.uniform([0.0, 0.0], [4.0, 3.0], (N, 2)) vel np.zeros((N, 2)) radius np.full(N, 0.03) for frame in range(240): # 重力加速度 vel g * dt # 位置更新 pos vel * dt # 地面碰撞反弹 floor_y 0.0 below pos[:, 1] floor_y radius pos[below, 1] 2 * (floor_y radius[below]) - pos[below, 1] vel[below, 1] -0.8 * vel[below, 1]注意看第14行到第16行的地面处理这里用了一个镜像反弹技巧求越界位置相对于边界的镜像位置而不是简单地把y坐标塞回边界值。这样能保证速度和位置的关系不产生新的能量突变。这段代码已经是完整的“粒子重力边界碰撞”物理系统用matplotlib.animation加个画图就能直接看到500个球在下落反弹。3.2 三种积分器怎么选数值积分的选择直接决定整个引擎会不会“爆”。我用过三种把各自特性说清楚显式欧拉是先更新位置再用新位置算加速度或者先算位置再更新速度两种变体在有弹簧这类力时会持续向系统注入能量模拟几十帧后整个场景就开始发疯。半隐式欧拉也叫symplectic Euler是先更新速度、再用新速度更新位置vel acc * dt pos vel * dt这个顺序的小改动带来质变因为对于自由落体和弹簧系统半隐式欧拉能保持系统能量上下波动而不是单调发散。这是游戏物理默认选择也是我所有后续代码的基础。Verlet积分用前一帧位置和当前帧位置做差分new_pos 2 * pos - old_pos acc * dt * dt它没有一个显式的速度量但天然适合约束很多的场景——布料、绳索、粒子堆叠。因为约束求解时只调整位置、从位置变化反推速度误差累积特性比直接操作速度好很多。选型建议初学阶段无脑用半隐式欧拉加固定时间步长。固定步长很关键可变步长会让物理系统在不同分辨率下表现不一致这是很多诡异抖动的根源。3.3 墙壁反弹与恢复系数边界处理除了镜像反弹还要引入恢复系数e它描述碰撞后法向速度保留多少。e1是完全弹性碰撞能量不损失e0是碰撞后法向速度清零物体贴住地面不动。实际场景一般取0.5到0.9之间代码里就是一个系数vel[below, 1] -0.8 * vel[below, 1]有个很容易被忽略的问题当球的落速极小时反弹速度也极小物体停不下来会一直在地面附近蠕动视觉上很烦。解决方法是给一个休眠阈值速度低于某个值直接清零threshold 0.05 low_speed np.abs(vel[below, 1]) threshold vel[low_speed, 1] 0.0这个细节在游戏手感里特别重要我后来做交互原型时几乎每个项目都用到。4. 刚体动力学与碰撞约束从“点”到“块”4.1 刚体的状态多出来的旋转自由度粒子只有两个自由度x和y刚体多一个旋转自由度θ。二维刚体的状态量变成五个位置(x, y)、速度(vx, vy)、角度θ、角速度ω。同样用数组组织pos np.zeros((N, 2)) vel np.zeros((N, 2)) angle np.zeros(N) # 物体朝向 omega np.zeros(N) # 角速度 I np.zeros(N) # 转动惯量转动惯量I是旋转动力学里的“质量”。实心圆盘对质心轴的转动惯量是0.5 * m * r^2它决定了同样力矩下角加速度有多大。刚体局部坐标系里的一个点r_local世界坐标计算公式是R np.array([[np.cos(angle[i]), -np.sin(angle[i])], [np.sin(angle[i]), np.cos(angle[i])]]) world_point pos[i] R r_local如果有一个力F作用在物体上但作用点不在质心它除了产生线加速度还会产生力矩τ。二维里的力矩是一个标量τ rx * Fy - ry * Fx其中(rx, ry)是从质心指向受力点的向量。这个公式在实现“推力器”“风帆”这类玩法时非常有用。4.2 碰撞检测Broad Phase 与 Narrow Phase碰撞检测是物理引擎最容易拖慢性能的部分。朴素做法是遍历所有物体两两比较复杂度O(N^2)500个球每帧要做12.5万次距离判断5000个球就是1250万次完全不可持续。工业级做法分成两阶段。Broad Phase快速排除明显不相交的物体对把候选对数量降到最小。最简单实用的是空间网格格子哈希把空间划分成边长略大于最大物体尺寸的均匀格子每个物体登记到自己所在的格子然后只检查同格子和相邻九格内的物体。因为物体密度通常不高平均复杂度能降到接近O(N)。Narrow Phase再精确判断候选对是否真的相交。球和球用距离平方和比较凸多边形用分离轴定理SAT——依次检查每条边的法线方向如果存在任何一个方向上两个图形的投影没有重叠就说明图形分离。这个定理是不相交判定逻辑很干净所有凸多边形碰撞的基础都是它。我在实际项目里的建议是先用朴素O(N^2)把逻辑跑通再用网格优化。过早优化会让代码复杂度翻倍掩盖了物理逻辑本身的问题。4.3 约束求解冲量、位置修正和一招稳定术碰撞发生后要做两件事阻止互相穿透、改变运动方向。经典做法是冲量法。两个球碰撞时沿法线方向的相对接近速度为vn碰撞后要变成方向相反、大小为e倍。解这一个方程得到法向冲量的公式j -(1 e) * vn / (inv_mass_i inv_mass_j)然后把冲量施加到两个物体上vel[i] j * inv_mass_i * n vel[j] - j * inv_mass_j * n这个公式在等质量弹性碰撞下正好退化成“交换法向速度”物理直觉完全对得上。但光改速度还不够离散时间步会让两个球在检测出碰撞时已经嵌入了一小段距离。只看速度不修正位置时间一长物体就会“沉”到地面里。所以我加了位置修正把重叠量按质量比例分摊给两个物体各推回一半质量大的动得少质量小的动得多。交给读者一个小技巧约束解算不要只跑一遍。很多引擎同一帧内会迭代2到4次约束修正每次把残留误差减半。迭代3次的效果和迭代20次接近而性能开销小得多。这也是为什么很多“堆叠箱子”演示里底部箱子看起来稳如泰山——其实是迭代约束的功劳。再进一步就是位置动力学PBD先根据力预测一个新位置然后把所有约束当成投影操作直接在位置上修正最后用修正后的位置反推速度。这套思路写布料、绳索、果冻这类软体模拟几乎是无脑好用的。4.4 提高模拟稳定性的几个工程细节经验的坑踩多了之后我习惯性地在刚体系统里加这些保护最大速度限制。每帧更新完速度后把所有速度幅值钳制在一个上限内比如20 m/s。这能防止离散碰撞检测漏检高速物体也能阻止数值溢出。阻尼项。每帧乘一个略小于1的系数比如vel * 0.999。它模拟空气阻力同时把数值误差缓慢压掉。休眠机制。速度低于阈值且连续数帧不变就把物体标记为休眠跳过它的积分和碰撞计算。大量堆叠刚体稳定后休眠能省下90%以上的计算量。物体的摩擦力近似。在碰撞速度修正里把切向速度乘以一个摩擦系数再保留一部分能产生粗糙表面的效果不需要引入复杂的摩擦锥模型。这些在引擎代码里往往只有几行但对结果的影响是决定性的。5. 性能优化向量化、预分配与Numba5.1 消灭for循环的两条路线Python物理引擎的性能瓶颈九成在for循环。消灭它们有两条路线。第一条是批量向量化——把外层循环去掉一次对整批数据做运算。比如施加外力代码从for i in range(N): force[i] ...变成force ...整个数组一次性算。或者把多个物体的运动统一写成矩阵运算一次执行。第二条是**“半向量化”**外层保留一个轻量循环循环内全部用向量操作。以球-球碰撞检测为例对于每个球i一次计算它和其他所有球的距离矩阵d pos - pos[i] # 广播所有球相对球i的偏移 dist_sq np.einsum(ij,ij-i, d, d) # 每个球到球i的距离平方 min_d radius[i] radius collide (dist_sq 0) (dist_sq min_d * min_d) idx np.flatnonzero(collide) # 找到真正碰撞的球索引这样外层循环N次内层全是C级别速度的向量操作比纯Python的两重循环快一到两个数量级而且可读性依然很好。5.2 减少内存分配的三个习惯NumPy数组加法a b会生成一个全新的临时数组物理引擎每帧要做大量这样的操作临时数组多了内存分配和释放本身就会拖慢速度。三个习惯能明显改善第一能原地操作就原地操作。vel g * dt这种是原地操作不产生新数组但如果写vel vel g * dt就产生新对象。第二带out参数计算。np.add(a, b, outc)把结果写入预分配的c省掉中间临时量。第三避免在循环里拼接数组。反复调用np.concatenate会让内存不断搬迁正确做法是先分配最大容量再对有效部分切片赋值。这三个习惯都改完之后我之前一个波形仿真项目性能提升了30%左右代码逻辑完全没变。5.3 Numba在不改物理代码的前提下再快一个数量级如果需要更大提升我推荐Numba——一个直接把Python函数编译成机器码的JIT库。它对NumPy的支持相当好很多Python循环代码只要加个装饰器就能提速几十倍。把积分和碰撞检测写成一个函数加上njitfrom numba import njit njit(cacheTrue) def physics_step(pos, vel, radius, inv_mass, dt): N pos.shape[0] vel g * dt # 球与球碰撞的纯Python循环 for i in range(N): for j in range(i 1, N): dx pos[i, 0] - pos[j, 0] dy pos[i, 1] - pos[j, 1] dist_sq dx * dx dy * dy min_d radius[i] radius[j] if dist_sq min_d * min_d and dist_sq 0: # 速度冲量计算 ... pos vel * dt这里有意思的地方在于Numba允许你写回“朴素的二重循环”因为编译后根本不再经过Python解释器。这比写一堆精妙的NumPy向量化技巧还快并且代码更直白。我自己处理几千个刚体碰撞时纯Python版跑2000帧要十几秒改成Numba后两秒以内稳跑。5.4 用一张表看三种实现方式的差距拿“600个小球堆积在箱子里模拟1000帧”做了一个很粗糙的基准三类写法的差距大约是实现方式骨架难度1000帧耗时我的旧笔记本可读性纯Python二重循环低22秒左右最好NumPy半向量化中3.5秒左右良好Numba JIT编译中0.8秒左右良好优化后的C高0.3秒左右较差Numba方案在代码量接近Python的前提下性能已经摸到了C的脚后跟。对于大多数教学、机器人训练和VR原型场景这个速度足够了。这套优化路径是我个人最喜欢的方向先追求逻辑清晰再考虑向量化最后上JIT每一步都有明确的收益。6. 实战500个球的可交互物理沙盒6.1 场景设计与参数设定理论讲完上一段真正能跑的代码。我设计了一个2D物理沙盒500个半径不等的球落入一个封闭矩形空间受重力影响彼此碰撞、和墙壁碰撞并反弹。关键参数取这些值时间步长1/120秒这样每帧画面有约8个物理子步视觉平滑且较稳定恢复系数0.85模拟出略带弹性的碰撞世界尺寸4米乘3米配合v0.05左右的球半径一个画面里能塞下足够的物体量。为什么是120Hz而不是60Hz因为物理系统对时间步长的敏感度和画面刷新率不是一回事。物理子步越小离散检测穿透的风险越低约束的收敛性越好。游戏引擎里把物理频率定在120Hz或240Hz很常见哪怕画面只刷新60次。6.2 完整代码实现我把它整理成一个完整可运行的脚本import numpy as np # 参数 N 500 DT 1.0 / 120.0 W, H 4.0, 3.0 G np.array([0.0, -9.8]) E 0.85 DAMP 0.999 # 状态 pos np.random.uniform([0.1, 0.1], [W - 0.1, H - 0.1], (N, 2)).astype(np.float64) vel np.zeros((N, 2)) radius np.random.uniform(0.03, 0.06, N) inv_mass np.ones(N) def step(): global pos, vel # 重力 空气阻尼 vel G * DT vel * DAMP # 球与球碰撞 for i in range(N): d pos - pos[i] # (N, 2) dist_sq np.einsum(ij,ij-i, d, d) min_d radius[i] radius collide (dist_sq 0) (dist_sq min_d * min_d) idx np.flatnonzero(collide) if idx.size 0: continue dist np.sqrt(dist_sq[idx]) n_ij d[idx] / dist[:, None] # 从i指向j的单位向量 rel_v vel[i] - vel[idx] vn np.sum(rel_v * n_ij, axis1) approaching vn 0 # 正在接近才处理 idx_app idx[approaching] if idx_app.size 0: continue n_app n_ij[approaching] vn_app vn[approaching] # 冲量 denom inv_mass[i] inv_mass[idx_app] j -(1 E) * vn_app / denom vel[i] np.sum((j * inv_mass[i])[:, None] * n_app, axis0) vel[idx_app] - (j * inv_mass[idx_app])[:, None] * n_app # 位置修正按质量比分摊重叠 overlap (radius[i] radius[idx_app]) - dist[approaching] w_i inv_mass[i] w_j inv_mass[idx_app] total_inv w_i w_j if total_inv 0: pos[i] np.sum((-overlap * w_i / total_inv)[:, None] * n_app, axis0) pos[idx_app] (overlap * w_j / total_inv)[:, None] * n_app # 推进位置 pos vel * DT # 墙壁约束与反弹 left pos[:, 0] radius right pos[:, 0] W - radius bottom pos[:, 1] radius top pos[:, 1] H - radius pos[left, 0] 2 * radius[left] - pos[left, 0] vel[left, 0] -E * vel[left, 0] pos[right, 0] 2 * (W - radius[right]) - pos[right, 0] vel[right, 0] -E * vel[right, 0] pos[bottom, 1] 2 * radius[bottom] - pos[bottom, 1] vel[bottom, 1] -E * vel[bottom, 1] pos[top, 1] 2 * (H - radius[top]) - pos[top, 1] vel[top, 1] -E * vel[top, 1] # 低速休眠阈值 vel[np.abs(vel) 0.001] 0.0 # 跑1200帧 for frame in range(1200): step()在Notebook里配合matplotlib.animation每50帧画一张散点图就能看到几百个球翻滚落稳的全过程。6.3 关键代码逐段讲解碰撞处理是整段代码的难点我拆开讲几个重点单位向量与接近判断。d pos - pos[i]计算的是“从球i指向其他球”的向量除以距离后就得到每个碰撞对象指向当前球i的单位向量n_ij。rel_v vel[i] - vel[idx]是球i相对碰撞对象的速度把它投影到n_ij上得到法向相对速度vn。只有vn大于0才说明两个球正在接近小于0说明正在分离无需处理——这个判断省掉了大量无用的速度修正。冲量施加。j -(1 E) * vn_app / denom计算出的是一个标量冲量沿着法线方向。球i收到的速度变化是j乘以自身的逆质量再乘方向向量碰撞对象则相反减去同样的量。因为这里是批量计算用了[:, None]扩展维度配合广播。位置修正的权重分配。重叠量必须消除但让谁动得多、谁动得少按逆质量比例分配。静态物体inv_mass为0时total_inv可能为0需要判断否则会除零。这个判断虽然小但少了它一旦混入静态物体引擎会瞬间崩掉。墙壁镜像反弹。墙壁部分没有用循环直接对布尔掩码索引的数组做赋值这正是NumPy风格的体现。镜像公式在边界之外也能正确处理不会卡死。6.4 向游戏、机器人、VR场景扩展沙盒能跑之后扩展方向已经很明朗。游戏方向把多边形的碰撞检测补上SAT加入关节约束旋转铰链、滑橇就可以做平台跳跃的谜题原型或者模拟弹弓、荡绳这些小玩法。Unity和Godot的物理引擎当然比这个完善但自己写一遍才能理解引擎参数背后的含义调起手感来也更准。机器人控制方向把动力学部分替换成机械臂的刚体链加上关节力矩就能和强化学习环境对接。我用Gymnasium自定义环境时核心的step函数里就是一套和上面结构完全相同的NumPy物理更新。批量状态处理用3D数组还能并行开多个环境的训练采样效率很高。VR培训方向沙盒里的碰撞反馈逻辑可以直接转成手势交互原型。比如模拟工具碰触虚拟物体时的震动强度用碰撞冲量大小作为反馈源的参数用户体验会非常自然。7. 常见问题排查与避坑实录7.1 环境类问题的速查表我把开箱时碰到的高频问题整理成表方便以后直接翻问题可能原因解决办法No module named numpy当前解释器不是虚拟环境VSCode里切换解释器或重新激活.venv安装卡在backend dependenciespip版本太旧尝试源码构建升级pip或--only-binary :all:强制用预编译包numpy版本不匹配报错Python版本太新、numpy太旧pip install --upgrade numpy代码在A电脑正常B电脑报错依赖版本不一致用requirements.txt锁定版本号运行中发现结果不确定用了随机数却没固定种子初始化时np.random.seed(42)7.2 物理模拟数值不稳定的几个典型症状症状一堆叠物体持续抖动、缓慢漂移。八成是约束迭代次数不够或时间步长太大。先把物理步长降到120Hz再把碰撞位置修正迭代次数提到4你会发现整个场景像被按住了。症状二快速飞行的球直接穿透墙壁。这是离散碰撞检测的经典缺陷——物体在一个时间步内从墙壁一侧运动到另一侧检测完全没碰到。解法有收紧时间步、给物体加连续碰撞检测、或者限制最大速度。我常用后两者组合。症状三整体能量越变越大物体会自己飘起来。这是积分器不稳定加碰撞处理给系统注入能量。检查是不是用了纯显式欧拉改成半隐式检查恢复系数是不是设成了大于1最后加上阻尼项兜底。症状四刚体旋转起来越来越快最后失控。通常是转动惯量算错或者力矩计算时把作用点向量用反了。二维里力矩τrxFy-ryFx方向正负搞反会让角速度往错的方向累加。7.3 我踩过的最贵的几个坑第一个坑是用逆质量表格处理约束时忘掉静态物体。第一次让箱子堆到地面后一个除零直接把整个进程送走。后来养成了习惯任何涉及inv_mass的除法先判断分母。第二个坑是为了追求全向量化把碰撞对全部构造成pair矩阵。1000个球会产生50万行pair中间矩阵瞬间吃掉几个G内存直接被系统警告。这是个很典型的“向量化思维过猛”案例真实场景应该用broad phase配合半向量化而不是把NPair全铺开。第三个坑是把物理步长和渲染帧率绑定。早期代码直接每帧调一次step结果在低帧率机器上物理明显变慢在高帧率机器上又太快。后来统一改成固定物理步长加累计器把渲染和物理完全解耦所有屏幕上的表现才一致起来。第四个坑比较深刻把所有物体集中到少量大数组之后调试几何问题变得极其痛苦。某个球的“位置”在数组里只是其中一行出了诡异行为没法直观看到。后来我自己习惯是核心计算用数组但排查问题时写一个def debug_obj(i): print(pos[i], vel[i], radius[i])之类的辅助函数或者临时只保留两三个物体跑最小复现场景。大多数情况下把问题缩小到最少重现条件原因一眼就能看出来。我个人在这些项目里的体会是物理引擎最有意思的不是把公式写对而是在“写对公式之后”的那个晚上——你盯着屏幕上几百个球自己推来推去发现它们有了沉甸甸的质感。PythonNumPy这个组合恰好能让你在最短时间内走到那个晚上。先把这套基础跑通后面不管是加关节、加流体还是加软体都是顺着同一套数组思维延伸出去的事。关于版本兼容性的问题建议换新机器时先跑一遍沙盒脚本确认环境一致再继续开发能省下一整天的排查时间。