我一直觉得很多人在听到“Python 写物理引擎”时第一反应是“这台机器怕不是要烧起来”。毕竟 Python 是出了名的解释执行、循环慢、GIL 锁死拿它去跑逐帧的物理模拟听起来确实像是在做性能自杀。但如果你从一开始就把所有物理状态组织成 NumPy 数组让热循环下沉到 C 层情况就完全不同了。这篇文章我想聊的就是如何用 Python NumPy 从零搭一个能用于游戏开发、机器人控制与 VR 交互原型的物理模拟引擎包括最核心的数值积分、碰撞处理、约束求解、性能优化以及一条条踩出来的经验。先说结论用 NumPy 写物理引擎本质不是“用 Python 写引擎”而是“用 Python 脚本当指挥中心让编译好的 C 代码干重活”。位置、速度、受力全部存在连续内存里一次加法直接操作整个数组手感上像在写 Python吞吐量却接近 C。这篇文章适合三种人想快速验证物理机制的游戏开发者、需要搭建简化动力学仿真环境的机器人研究者、以及所有好奇“NumPy 为什么能加速”的人。下面我从最底层的设计思路开始逐步展开。1. 为什么物理引擎要拥抱 NumPy一场关于“慢语言”的误会1.1 “慢语言”与“快数组”被误解的事实Python 慢这个结论本身没错。问题在于很多人的慢是慢在写法上——三层 for 循环嵌套逐个粒子更新状态。每次循环都要做解释器调度、对象属性查找、数值装箱拆箱跑 100 万次再快的机器也扛不住。但 NumPy 的核心逻辑是预编译的 C 和 Fortran 代码它把“对每个元素执行同一操作”的任务批量完成。用 NumPy 写物理引擎时你的循环通常只发生在以下两处一是每个时间步的全局更新二是少量约束迭代。真正逐粒子的操作全部通过向量化表达式完成。我用一个简单例子验证过。100 万个粒子每个粒子有位置、速度、外力做一次最基础的“受力更新速度速度更新位置”import numpy as np import time n 1_000_000 positions np.random.rand(n, 3) velocities np.random.rand(n, 3) * 0.1 forces np.random.rand(n, 3) mass np.ones(n) dt 0.016 # 纯 Python 写法先转列表再逐粒子更新 t0 time.perf_counter() pos_list positions.tolist() vel_list velocities.tolist() f_list forces.tolist() for i in range(n): ax f_list[i][0] / mass[i] ay f_list[i][1] / mass[i] az f_list[i][2] / mass[i] vel_list[i][0] ax * dt vel_list[i][1] ay * dt vel_list[i][2] az * dt pos_list[i][0] vel_list[i][0] * dt pos_list[i][1] vel_list[i][1] * dt pos_list[i][2] vel_list[i][2] * dt t1 time.perf_counter() print(纯 Python 耗时:, t1 - t0, 秒) # NumPy 写法 t0 time.perf_counter() acceleration forces / mass[:, np.newaxis] velocities acceleration * dt positions velocities * dt t1 time.perf_counter() print(NumPy 耗时:, t1 - t0, 秒)我自己机器上测出来的差距大概是几十倍到上百倍。这不是说 Python 开发者笨而是“命令式循环”和“数组表达式”属于完全不同的执行模型。理解这一点是设计高性能物理模拟的第一个门槛。1.2 物理世界天然是列式数据很多初学者写物理模拟第一直觉是建一个 Particle 类里面放 x、y、z、vx、vy、vz然后 new 一万个对象塞进列表。这在 Python 里是灾难。每个对象都有独立的 dict 存储属性内存碎片化严重遍历列表等于不停地在指针之间跳来跳去CPU 缓存利用率极低。NumPy 的思路是“列式存储”不再维护一背包粒子对象而是维护几个大数组。x 是全体粒子 x 坐标的数组y 是全体粒子 y 坐标的数组或者更常见的直接用一个(N, 3)的二维数组存所有位置。数组里每个元素紧挨着排布在连续内存上做一次加法时CPU 可以顺序读、顺序写缓存命中率完全不同。用生活化类比就是你有一百个快递要送一个地址一张纸挨个翻查效率低不如拉一张大表所有门牌号都在同一列里一列一列处理。物理模拟本质上是批量计算天然适合这种“拉表”式数据组织。1.3 向量化带来的第一道性能红利向量化还意味着你可以直接利用 CPU 的 SIMD 指令。现代处理器支持单条指令同时处理多个浮点数NumPy 底层会在合适的时机做这种优化。你不需要手动写汇编只需要把你的计算表达成“数组加数组”“数组乘标量”这种形式剩下的交给底层。当然向量化不是银弹。它适合“结构一致、计算同构”的任务而不适合“大量分支、稀疏操作”的任务。物理模拟正好是前者的典型代表每帧每个粒子都执行同样的牛顿第二定律无非是力不同、质量不同完全可以用同一条数组表达式覆盖。2. 状态量的数组化设计把物理场装进内存2.1 最小状态集位置、速度与力一个刚体或粒子系统核心状态量其实就三个位置、速度、外加合力。质量可以视为常数。用 NumPy 组织起来极其简单N 10000 positions np.zeros((N, 3), dtypenp.float64) velocities np.zeros((N, 3), dtypenp.float64) forces np.zeros((N, 3), dtypenp.float64) mass np.full(N, 1.0, dtypenp.float64)这里有几个细节值得注意dtype默认就是float64显式写出来是为了提醒自己“一切都是连续内存块”。(N, 3)的排布有一个额外好处调用forces[:, 1] ...可以直接处理所有粒子的 Y 分量比如施加统一重力。力的数组每帧都要清零重建所以forces[:] 0.0会比forces np.zeros_like(forces)更好因为前者原地清零不触发新内存分配。我见过不少新手在循环里反复执行positions np.array([...])每帧创建新数组导致内存申请和释放的 overhead 非常大。正确做法是初始化一次之后所有更新都靠和*完成保持数组身份不变。2.2 力模型不是循环是数组表达式最常见的三种力——重力、线性阻力、弹簧力——在 NumPy 里都只有几行。gravity np.array([0.0, -9.8, 0.0]) k_drag 0.1 forces[:] 0.0 forces[:, 1] mass * gravity[1] # 重力F m * g forces - k_drag * velocities # 线性阻力F -k * v弹簧力稍微复杂一点因为需要计算粒子对之间的相对位置。做法是维护两个索引数组src和dst分别表示弹簧两端的粒子编号然后通过高级索引一次性取出两端位置src np.array([0, 1, 2, ...], dtypenp.int32) dst np.array([1, 2, 3, ...], dtypenp.int32) delta positions[src] - positions[dst] dist np.linalg.norm(delta, axis1, keepdimsTrue) direction delta / (dist 1e-12) # 加一点保护避免除零 spring_force_magnitude -k * (dist - rest_length) forces[src] direction * spring_force_magnitude forces[dst] - direction * spring_force_magnitude这里每一步都是对整个弹簧集合起作用的即使有几万条弹簧也只是一次数组运算而不是一个 Python 循环。2.3 心智模型对数组说话而不是对粒子说话写熟之后你会形成一种新的心智模型不再问“这个粒子受到哪些力”而是问“这一批粒子的某个分量整体满足什么关系”。比如“所有低于地面的粒子把 Y 坐标设回 0”——这不是 if 语句是布尔掩码。“所有弹簧的当前长度小于自然长度就施加推力”——这不是分支是符号运算。“把所有粒子的速度同时加上重力加速度乘 dt”——这不是循环是一行乘法加法。这个转变对工程实现很重要。只要你开始像操作电子表格一样操作物理数组写物理引擎的难度会大幅下降。3. 积分器选型半隐式欧拉为什么是默认答案3.1 三种积分器的对比物理引擎的核心循环就是不断积分知道加速度算速度算位置。常见选择有三种积分器更新顺序稳定性适用场景显式欧拉先算位置再算速度都基于旧状态能量容易漂移弹簧系统会越跑越“兴奋”简单的学习演示不推荐用于引擎半隐式欧拉semi-implicit Euler先更新速度再用新速度更新位置相对稳定系统呈轻微耗散趋势游戏开发、机器人仿真的默认选项速度 Verlet先更新半速再更新位置再用新加速度更新半速对振荡系统更稳定能量守恒性更好分子动力学、布料模拟为什么显式欧拉不稳定因为它用“旧速度”推动“新位置”相当于每个时间步都在给系统强行注入误差在弹簧这种周期性运动里误差会累积成震荡发散。半隐式欧拉先用新速度推动位置等于把一部分未来信息提前纳入计算无形中给系统加了阻尼。3.2 半隐式欧拉的五行实现实现极简def integrate(positions, velocities, forces, inv_mass, dt): acceleration forces * inv_mass[:, np.newaxis] velocities acceleration * dt positions velocities * dt forces[:] 0.0注意这里用了inv_mass逆质量而不是直接forces / mass。原因后面会讲。这五行的顺序不能乱一定是先速度后位置。如果你把顺序写反就退化成显式欧拉。3.3 步长与刚度的平衡术积分步长dt不是随意定的。物理系统有一个特征频率比如一根弹簧刚度k100质量m1角频率ωsqrt(k/m)10 rad/s周期约 0.63 秒。为了保证数值稳定性dt至少要小于振荡周期的几十分之一经验上取1/(10ω)以下也就是 0.01 秒级别。弹簧越硬dt必须越小。这也是为什么布料模拟里弹簧刚度过大时布会直接炸开——不是物理规则错了是数值方法跟不上。我的经验法则是先按目标帧率选dt比如 60 FPS 则dt1/60≈0.0167然后反推系统里允许的最大弹簧刚度。如果刚度过高要么降低刚度要么减少dt在游戏原型里可以直接做“每帧内部跑多个物理子步”substeps 4 sub_dt dt / substeps for _ in range(substeps): integrate(...)这样物理精度提高了但计算量也翻倍。怎么取舍取决于你对实时性的要求。4. 碰撞检测与响应复杂度才是真正的坎4.1 为什么碰撞是必需品一个只会飘粒子的“物理引擎”本质上只是个动画播放器。真正的物理交互来自碰撞球落在地上弹起、布料搭在桌上、机械臂末端触碰物体。在游戏和机器人场景里碰撞决定了一切“接触感”。它也是从“数组计算”走向“算法设计”的关键分水岭。4.2 朴素最近邻与广播的威力最容易想到的碰撞检测是遍历所有粒子对判断距离是否小于半径之和。NumPy 里可以一行广播得到所有粒子对距离diff positions[:, np.newaxis, :] - positions[np.newaxis, :, :] dist_sq np.sum(diff * diff, axis2) np.fill_diagonal(dist_sq, np.inf) pairs np.argwhere(dist_sq collision_threshold_sq)这写法优雅、直观但只能在小规模下使用。原因无他复杂度是 O(N²)内存也是 O(N²)。当 N10000 时dist_sq是一个 10000×10000 的 float64 矩阵占用约 800 MBdiff更是要到 2.4 GB。N100000 时数据量直接到 TB 级内存立刻爆炸。4.3 空间哈希把复杂度拖回地面工程上更实用的方案是空间哈希。思路是把三维空间切成固定大小的网格每个粒子只和自己所在格子及相邻格子里的粒子做碰撞检测。这样平均复杂度降为 O(N×每个格子内粒子数)在粒子稀疏分布时接近线性。def build_spatial_hash(positions, cell_size): cell_coords np.floor(positions / cell_size).astype(np.int32) buckets {} for idx, cell in enumerate(cell_coords): key (cell[0], cell[1], cell[2]) buckets.setdefault(key, []).append(idx) return buckets检测碰撞时遍历每个格子再检查自身和周围 26 个相邻格子的粒子对。这个 Python 循环仍然存在但每个格子里粒子数很少循环总量可控。如果还想继续提速可以把这个函数用 Numba 的njit编译后面会提到。4.4 碰撞响应修位置、反弹速度找到碰撞对之后关键是处理“穿透”。大多数情况下模型已经发生了微小穿透只靠速度反弹是不够的。两步走for i, j in contact_pairs: delta positions[j] - positions[i] dist np.linalg.norm(delta) if dist 2 * radius: normal delta / (dist 1e-12) overlap 2 * radius - dist positions[i] - normal * overlap / 2 positions[j] normal * overlap / 2 rel_vel velocities[j] - velocities[i] vn np.dot(rel_vel, normal) if vn 0: impulse -vn * (1 restitution) velocities[i] - normal * impulse / 2 velocities[j] normal * impulse / 2先修正位置再用相对速度沿法线方向做弹性反弹。restitution是恢复系数0 表示完全非弹性1 表示完全弹性。工程里通常取值 0.2~0.8既能表现碰撞又不至于让物体弹个没完。这个循环是 Python 层逐对处理的对几百对碰撞来说没问题。当碰撞对上千时就要考虑向量化或 Numba。物理引擎的优化通常从这里开始。5. 性能优化实操从“能跑”到“跑得快”5.1 消灭 Python 层循环性能优化第一条检查所有模拟热路径里的 Pythonfor能换成数组表达式的坚决换掉。最典型的例子是边界碰撞处理。慢写法for i in range(N): if positions[i, 1] 0: positions[i, 1] 0 velocities[i, 1] * -0.5向量化写法below_ground positions[:, 1] 0 positions[below_ground, 1] 0 velocities[below_ground, 1] * -0.5前者是逐粒子判断后者是一次布尔掩码加两次批量赋值。效果完全一样但后者在处理百万粒子时优势巨大。条件分支很多时候也可以用np.where或乘法遮挡表达# 阻尼力只在速度超过阈值时施加 damping np.where(np.abs(velocities) threshold, 0.0, k_drag) forces - damping * velocities5.2 预分配内存一次申请反复使用Python 的 GC 会对小对象的创建销毁造成额外开销。哪怕 NumPy 数组底层不参与 GC但每帧np.zeros_like新开一个数组依然要付出系统内存分配的成本。正确姿势是所有中间量在初始化阶段就申请好后续用np.multiply、np.add的out参数原地写入。delta_v np.empty_like(velocities) np.multiply(acceleration, dt, outdelta_v) velocities delta_v同样地forces[:] 0.0与forces np.zeros_like(forces)的区别就在于此前者复用原内存后者重新分配。5.3 数据布局AoS 还是 SoA我在前面已经提到了列式存储。这里展开说下两个方案AoSArray of Structuresparticles np.zeros((N, 3))每个粒子的三个分量挨在一起。写起来直观提取单个粒子的速度容易。SoAStructure of Arraysvx np.zeros(N); vy np.zeros(N); vz np.zeros(N)每个分量单独一个数组。物理模拟里我倾向于 SoA 或“类 SoA”的二维数组因为向量化运算时对连续分量的批量操作更容易命中缓存。如果你只关心位置和速度的数组表达式(N, 3)已经够用。但当粒子数量极大、并且你频繁只操作某一个分量时拆开存会有额外收益。大家可以根据需要做实验不必盲从一个方案。我自己在 50 万粒子规模下测过SoA 的带宽利用率大约比 AoS 高 20%-30%但代码可读性差一些。原型阶段用(N, 3)性能优化阶段再根据热点迁移也不迟。5.4 当 NumPy 不够用Numba 与 JAX 的接力如果算法里确实存在需要大量分支的循环比如碰撞对处理、约束迭代NumPy 很难优雅表达代码反而更乱。这时候我用 Numba。它的njit装饰器可以直接把 Python 函数编译成机器码对循环的加速可以达到几十倍from numba import njit njit def satisfy_constraints_nb(positions, constraint_pairs, rest_lengths, iterations): for _ in range(iterations): for idx in range(constraint_pairs.shape[0]): a, b constraint_pairs[idx] delta positions[b] - positions[a] dist np.sqrt(delta[0] ** 2 delta[1] ** 2 delta[2] ** 2) if dist 1e-12: continue correction (dist - rest_lengths[idx]) / dist positions[a] delta * correction * 0.5 positions[b] - delta * correction * 0.5写法看起来像 Python实际跑起来接近 C。另一个选择是 JAX它的vmap、jit和 GPU 支持在科学计算里很有潜力但引入的编译依赖和心智成本更高。我的建议是先用 NumPy 把模型搭通再针对热点函数上 Numba不要一开始就上重型工具。6. 实战一个可运行的粒子布料模拟引擎6.1 布料建模布料可以抽象成一个网格粒子系统每个网格点是粒子相邻粒子之间用弹簧连接。为了简单又不失去布料的基本特性只需要两类弹簧水平方向连接同行相邻点垂直方向连接同列相邻点。再加固定点布的悬挂点和重力就是一个最早期的布料模拟雏形。6.2 约束迭代与弹簧力的区别很多人会用上一章提到的弹簧力来做但布料模拟里更推荐直接用“距离约束 迭代求解”。区别在于弹簧力是“柔”的刚度不够时布会像橡皮泥刚度过高时数值爆炸距离约束是“硬”的每次直接把两个粒子的距离拉回自然长度并且通过多轮迭代让所有约束相互协调。实现方式是对每条弹簧算出当前两端距离和自然长度的偏差按两端质量的倒数比例把位置修正掉。固定点质量设为 0不参与移动。整个算法的核心就是“位置投影”。6.3 完整代码与调参记录下面是一个 20×20 粒子布料的完整实现跑起来可以看到垂坠效果import numpy as np ROWS, COLS 20, 20 DT 0.005 SUBSTEPS 4 ITERATIONS 5 grid_x, grid_z np.meshgrid(np.linspace(0, 1, COLS), np.linspace(0, 1, ROWS)) positions np.stack([grid_x.ravel(), np.zeros(ROWS * COLS), grid_z.ravel()], axis1).astype(np.float64) velocities np.zeros_like(positions) mass np.ones(ROWS * COLS) # 固定布料左上角和右上角两个悬挂点 mass[0] 0.0 mass[COLS - 1] 0.0 inv_mass np.zeros_like(mass) safe mass 0 inv_mass[safe] 1.0 / mass[safe] # 生成弹簧约束 constraint_pairs [] for r in range(ROWS): for c in range(COLS - 1): a r * COLS c b r * COLS c 1 constraint_pairs.append((a, b)) for r in range(ROWS - 1): for c in range(COLS): a r * COLS c b (r 1) * COLS c constraint_pairs.append((a, b)) constraint_pairs np.array(constraint_pairs, dtypenp.int32) rest_lengths np.linalg.norm( positions[constraint_pairs[:, 0]] - positions[constraint_pairs[:, 1]], axis1 ) def apply_forces(forces): forces[:] 0.0 forces[:, 1] mass * (-9.8) forces - 0.05 * velocities def satisfy_constraints(): for _ in range(ITERATIONS): delta positions[constraint_pairs[:, 0]] - positions[constraint_pairs[:, 1]] dist np.linalg.norm(delta, axis1) correction (dist - rest_lengths) / (dist 1e-12) wa inv_mass[constraint_pairs[:, 0]] wb inv_mass[constraint_pairs[:, 1]] total wa wb factor_a correction * wa / (total 1e-12) factor_b correction * wb / (total 1e-12) positions[constraint_pairs[:, 0]] - (delta.T * factor_a).T positions[constraint_pairs[:, 1]] (delta.T * factor_b).T for step in range(300): sub_dt DT / SUBSTEPS for _ in range(SUBSTEPS): forces np.zeros_like(positions) apply_forces(forces) acceleration forces * inv_mass[:, np.newaxis] velocities acceleration * sub_dt positions velocities * sub_dt satisfy_constraints() # 简单地面碰撞 below positions[:, 1] 0.0 positions[below, 1] 0.0 velocities[below, 1] * -0.3代码里最值得注意的是(delta.T * factor_a).Tdelta是 (K, 3)factor_a是 (K,)要按粒子的行分别缩放需要先转置再乘再转置回来。当然用delta * factor_a[:, np.newaxis]更直观我这里展示的是另一种写法两者等价。调参经验ITERATIONS越大布越“刚”5 次能满足多数视觉需求SUBSTEPS越大系统越稳定但每帧耗时成倍增加DT0.005时布料下坠自然如果调大突然飞起来优先看是不是固定点没设对其次降低DT。7. 从原型到产品游戏、机器人和 VR 里的实际用法7.1 游戏开发Python 原型与 C 落地的分工很多人会问游戏引擎内部都有现成的物理后端为什么还要自己写答案是“原型验证”。当你设计一个新的关卡机制比如“玩家拉动绳索导致吊桥翻转”只想快速验证手感没必要在 Unity 或 Unreal 里折腾可视化脚本、物理材质、骨骼约束。用 Python NumPy 跑一个简化模型几小时就能把参数范围和玩法机制摸清楚再回到 C 里按同样的逻辑实现。我自己做过一次绳索物理验证Python 版本跑了 5000 个约束点在交互式调节参数时依然能实时更新。这个速度对原型足够。7.2 机器人控制仿真先行真机微调机器人领域更看重的是“控制策略验证”。直接在真机上试算法风险大、成本高。常见的做法是搭一个简化的刚体/软体仿真环境把状态向量给控制器控制器输出力矩再反馈回仿真。Python NumPy 可以快速搭出这类闭环。PyBullet 这类成熟仿真器底层也是 C但你能用 NumPy 自己搭一个“极简仿真器”专门针对某个机械臂的特定关节动力学这在研究侧非常常见。7.3 VR 交互16 毫秒的苛刻预算VR 的刷新率普遍是 90Hz 或 120Hz意味着留给物理计算的时间往往不到几毫秒。纯 Python 引擎做重型刚体模拟不现实但做简单的物体交互抓取、投掷、震动反馈是可行的。关键在于物理渲染的物体数量要克制碰撞检测用空间哈希而不是广播dt要固定并用子步控制稳定性。我见过一个 Demo用 NumPy 驱动 200 个刚体的简单碰撞再配合 GPU 渲染在 VR 原型机里能稳定跑满 90 帧。7.4 我踩过的那些坑最后集中整理几个我在实际调试中踩得最深、也最容易被新手忽视的坑固定点质量设为 0但惯性计算时忘了处理。forces / mass直接除零得到 NaN然后整个物理世界瞬间变成“一锅粥”。正确做法是像本文代码里一样用inv_mass并先屏蔽质量为零的点。dt 不一致导致结果完全不可复现。游戏渲染帧率是波动的物理模拟必须用固定dt否则同样的场景跑两次结果不一样。控制逻辑也容易在帧率波动时失稳。朴素碰撞检测的内存爆炸。小规模用广播很爽规模上去了要用空间哈希不然代码还没跑到响应阶段就 OOM。约束迭代时除以距离容易除零。两个粒子完全重合时dist0修正量变成 NaN。所有涉及范数的除法都要加一个小量1e-12做保护。NumPy 版本与随机数种子。如果你在模拟里用了np.random不同 NumPy 版本的默认随机算法可能不同导致随机初始条件不一致。需要精确复现时明确指定np.random.default_rng(seed)。我个人的习惯是任何涉及物理模拟的实验代码第一步先写一个“烟雾测试”——放两个粒子、一根弹簧、一块地面让系统跑几百步看能量是否单调或不合理地增长。这种小测试能过滤掉 90% 的积分和约束实现的低级错误比直接在完整场景里排查省事得多。物理引擎这东西表面上是数学实际写起来全是工程细节。把数组组织好把算法复杂度降下来再根据热点逐步优化这条路走通之后你会发现在很多需要快速验证想法的场合Python NumPy 远比想象中能打。