写这篇文章之前，我先把话说在前面：市面上的物理引擎已经做得足够好了，无论是闭源的商业引擎、开源的物理库，还是游戏引擎内置的物理模块，都能满足大多数场景。但你如果真的想搞懂“碰撞检测为什么这么写”“堆叠为什么总会抖”“摩擦力为什么有时候像抹了油”，唯一的办法就是自己从零写一个足够简单的实现，把每一行代码和背后的物理公式对上。这篇帖子就带你走完这个过程，最终你会得到一个能跑、能看、能继续扩展的二维物理引擎核心。全程不依赖任何第三方图形库，只要有一个能跑 Python 的环境，代码就能复现。

先说清楚目标：我们要实现的不是 Box2D 那样的工业级物理引擎，恰恰相反，它足够“简易”——支持圆形和轴对齐矩形，支持重力、冲量碰撞响应、库仑摩擦、基础堆叠和穿透修正。别小看这点能力，市面上不少独立游戏和原型 Demo 的物理效果就是从这个量级起步的。读完并跑通代码后，你会对刚体动力学、碰撞流形、迭代求解这些概念有完全不一样的体感。

1. 为什么需要亲手写一个“简易物理引擎”：目标、边界与引擎全景

动手之前得先想明白一件事：你为什么要造这个轮子。我的理由是，物理引擎表面上是一个 API，本质上是一套数学算法。你在使用现成引擎时遇到的大部分诡异问题——物体抖、穿透、突然弹飞——都源于对底层算法缺乏直觉。自己写一遍，你会建立起一套“物理引擎是在逼近真实，而不是计算精确解”的判断标准，后面再去调参、优化、选型，思路会清晰得多。

1.1 这个引擎要解决什么问题

我们针对的场景非常具体：一个二维世界，里面有圆球、方块、一面地面墙，重力往下拉，物体掉下来会互相碰撞、会堆叠、会滑动。不需要流体、软体、关节、粒子，也不需要连续碰撞检测。把这些功能砍掉之后，核心就只剩三件事：

• 碰撞检测：判断两个形状有没有相交，以及相交的深度、法线、接触点。
• 动力学求解：把重力和其它力转化为速度和角速度的变化。
• 碰撞响应：用冲量把相互穿透的物体推开，同时让速度发生合理的跳变。

这三件事恰好对应了物理引擎中最经典的三层结构，也是很多引擎源码里模块划分的原型。从工程实现的角度看，一个“简易物理引擎”的骨架大致如下面这个表：

这不是唯一的划分方式，但作为一篇教程的骨架，它足够清晰。后面的每一章都会落到这张表上的具体位置。

1.2 明确“不做”的部分

写代码之前，还有必要把边界划清楚。我们不做旋转矩形（OBB）之间的通用分离轴检测，只做轴对齐矩形（AABB）；不做关节约束，不做马达，不做流体；固定时间步长而不是自适应步长；用冲量法而不是基于力的约束求解。

这些简化是刻意的。轴对齐矩形的碰撞检测实现起来远比旋转矩形简单，但核心思路——分离轴定理——是一样的。先跑通 AABB 版本，再去扩展到 OBB，路径要平滑得多。这就好比你不会一开始就给引擎加上全套优化，先把一辆能跑的车造出来，再谈悬挂和空气动力学。

把目标订到这么具体还有一个额外的好处：一旦 Demo 跑通，你能亲眼看到“为什么堆叠需要迭代求解”“为什么穿透修正不能直接改位置”，因为你会亲手写出那版“会塌的代码”和“修好的代码”，这种对比比任何理论文章都来得直接。

2. 地基工程：矢量的加减乘除与刚体状态表示

物理引擎里 80% 的运算发生在二维平面上：力的合成、速度的分解、方向的判定。所以第一个要写好的部分就是一个可靠的二维矢量类。它不需要多高性能，但要语义清晰、不犯错。我在自己的实现里一般只写六七个方法，因为它们几乎构成了后面所有几何算法的最小集合。

2.1 一个够用的 Vec2 实现

下面这段代码是完整的二维向量类，我会把它作为整个引擎的数学地基。你注意看cross方法，在二维里它返回的是一个标量，这个标量的大小等于 |a||b|sinθ，符号表示旋转方向。这个量在冲量计算里非常关键，它会反复出现。

import math class Vec2: __slots__ = ("x", "y") def __init__(self, x=0.0, y=0.0): self.x = float(x) self.y = float(y) def __add__(self, other): return Vec2(self.x + other.x, self.y + other.y) def __sub__(self, other): return Vec2(self.x - other.x, self.y - other.y) def __mul__(self, scalar): return Vec2(self.x * scalar, self.y * scalar) __rmul__ = __mul__ def dot(self, other): return self.x * other.x + self.y * other.y def cross(self, other): return self.x * other.y - self.y * other.x def length_sq(self): return self.x * self.x + self.y * self.y def length(self): return math.sqrt(self.length_sq()) def normalized(self): length = self.length() if length < 1e-12: return Vec2(1.0, 0.0) inv = 1.0 / length return Vec2(self.x * inv, self.y * inv) def perpendicular(self): # 逆时针旋转90度：(-y, x) return Vec2(-self.y, self.x)

为什么这么设计？核心原因有三个。第一，所有方法都返回新对象，避免别名修改导致的隐蔽 bug；第二，length_sq专门用于距离比较，避免无谓的开平方，这在碰撞检测里是高频操作；第三，perpendicular方法后面摩擦和切向计算都要用它，单独抽出来比每次都现写Vec2(-v.y, v.x)要清晰得多。

2.2 刚体的状态字段：为什么要存质量倒数

有了矢量，下一步是定义刚体。初学者最容易忽略的一个设计是：不要直接存mass，要存invMass（质量的倒数）。原因很实际：物理引擎里大量公式形如“速度变化量 = 冲量 × invMass”，单位质量的倒数能直接把除法变成乘法，效率高一点只是其次，更重要的是静态物体（质量无限大）可以统一表示为invMass = 0，代码里不需要到处写if static分支。

同理，转动惯量的倒数invI也这样处理。转动惯量本身的计算随着形状不同而不同，圆是0.5 * m * r²，矩形是m * (w² + h²) / 12。但在简易引擎里，我建议先让圆的invI = 0，也就是不参与旋转，这样冲量公式可以先跑通最简版本，后面要开放旋转时再把对应项加回来。

刚体的完整状态我用这样一组字段：

class Body: def __init__(self, shape, position, mass=1.0, restitution=0.3): self.shape = shape # "circle" 或 "aabb" self.shape_data = ... # 圆的半径或矩形的半宽半高 self.position = position # 质心位置 self.velocity = Vec2() # 线速度 self.angle = 0.0 # 当前角度（AABB模式下恒为0） self.angular_vel = 0.0 # 角速度 self.force_acc = Vec2() # 每帧重置的力累加器 self.inv_mass = 0.0 if mass <= 0.0 else 1.0 / mass self.inv_i = 0.0 # 本例中默认不旋转 self.restitution = restitution # 恢复系数

注意force_acc的使用习惯：物理引擎里力通常是每帧清零后重新累加的，这让你可以在一个帧内给物体加多个力，然后在积分阶段一次性应用。如果不这么做，上一帧的力会残留下来，导致物体越跑越奇怪。

2.3 一个容易被忽略的约定：所有速度都必须相对于质心

刚体的线速度是相对质心的，角速度是绕质心的。在很多入门实现里，这两者是分开处理的，但如果你的接触点不在质心位置，速度场其实应该是：

接触点的实际速度 = 质心速度 + 角速度 × (接触点 - 质心) 的垂直向量

这句话现在不用完全理解，只要记住一点：冲量计算时，如果你忽略了接触点到质心的偏移，那么你算出来的相对速度就是不准确的，物体会出现“明明撞上了却不转”的诡异表现。我们 Demo 里的圆形物体invI = 0可以暂时跳过这个计算，但概念上要留个印象。

3. 运动积分：半隐式欧拉为什么是默认选项

物理引擎的每一帧要回答一个问题：给定当前速度、加速度和时间步长，物体下一秒该在哪里。这个过程的学名叫“数值积分”。你可能听说过 RK4、Verlet、各种高阶积分器的名字，但在绝大多数游戏物理引擎里，实际用的是最朴素的一阶方法——半隐式欧拉。

3.1 显式欧拉与半隐式欧拉的差距

显式欧拉是最直觉的写法：

位置 += 速度 * dt 速度 += 加速度 * dt

它先更新位置，再用旧速度叠加。问题是这种顺序会让系统能量漂移，典型表现是弹簧越蹦越高、球越弹越远，因为位置更新时用的速度已经“过时”了。半隐式欧拉把顺序反过来：

速度 += 加速度 * dt 位置 += 速度 * dt

一个顺序的调换，带来的稳定性提升是质的飞跃。位置更新用的是刚刚计算过的新速度，相当于显式地利用了未来的信息，能量漂移要小得多。你可以拿一个最简单的弹簧系统做实验，分别用两种顺序跑 10 秒，显式欧拉早就炸了，半隐式欧拉还在规规矩矩地振荡。这也是为什么工业级引擎在默认积分器上惊人一致地选择了半隐式欧拉。

简易引擎的定位决定了半隐式欧拉就是正确答案。实现同样简单得惊人：

def step(body, gravity, dt): # 力 -> 加速度 -> 速度（半隐式第一步） acceleration = body.force_acc * body.inv_mass body.velocity += acceleration * dt body.velocity += gravity * dt # 速度 -> 位置（半隐式第二步） body.position += body.velocity * dt body.force_acc = Vec2() # 清零力累加器

其中gravity我一般设为Vec2(0, -9.8)，朝下为负 y 方向。这样一种简单的推进方式，可以看作整个物理世界的“时钟”。

3.2 固定时间步长：为什么不能直接每帧调一次

很多初学者做引擎时会在渲染循环里直接调用step(dt)，而dt是上一帧的实际耗时。这个做法有个隐蔽问题：dt是波动的，而物理模拟对时间步长非常敏感。同一个场景，60 FPS 和 30 FPS 跑出来的结果可能完全不同，因为碰撞求解对步长的依赖是线性的。

正确做法是用一个累加器把渲染帧率与物理步长解耦：

dt_fixed = 1.0 / 60.0 accumulator = 0.0 def frame_update(frame_dt): global accumulator accumulator += min(frame_dt, 0.25) # 防止死亡循环 while accumulator >= dt_fixed: physics_step(dt_fixed) accumulator -= dt_fixed

每次physics_step都用固定的dt_fixed，模拟结果才可复现。这个技巧看起来简单，但它直接影响“堆叠会不会抖、碰撞会不会穿透”的体验。建议从第一天就把这个结构搭好。

3.3 积分顺序的细节

一个完整物理帧的积分顺序，我按下面这个顺序执行：

1. 合并所有外力到force_acc；
2. 用半隐式欧拉更新所有物体的速度；
3. 更新所有物体的位置和角度；
4. 碰撞检测，生成接触对列表；
5. 迭代求解所有接触对的冲量、摩擦；
6. 位置修正，把穿透拉回阈值内。

很多人会把碰撞响应放在速度更新之前，结果撞完还会穿透。顺序问题的本质在于：碰撞响应处理的是“最新位置上的速度关系”，所以检测应当基于最新的位置，求解应当作用于最新的速度。这篇博文后面所有代码都遵循这个顺序，建议你也固定下来，别轻易换。

4. 碰撞检测：从 Broad Phase 到 SAT 分离轴原理

碰撞检测是物理引擎里工作量最大、也最容易出错的部分。它分了两个阶段：Broad Phase 是粗筛，目标是快速排除不可能相交的物体对；Narrow Phase 是精算，对候选物体对做精确的几何测试，输出接触流形。

4.1 Broad Phase：先花小代价排除大多数

如果场景里的物体数量少于三四十个，一种足够简单的 Broad Phase 是直接暴力遍历所有物体对，每对先用包围盒做一次快速排斥测试。即使这样，也能过滤掉大量不相邻的物体。

包围盒用 AABB 表达就够了，每个刚体维护一个min_x, max_x, min_y, max_y。两个物体要相交，必须满足各自 AABB 在 x 和 y 方向都重叠。这里的优化点是：AABB 可以随着物体移动每帧重算，也可以惰性更新。简易引擎直接每帧重算即可，因为刚体数量不多。

如果想更激进一点，可以用一个均匀网格或者空间哈希表把世界划分成格子，每帧只检测同格子和相邻格子里的物体对。这个小优化大概几十行代码，能将 n² 的复杂度降为近似 n。对于二维场景尤其有效。最怕的就是物体都挤在一个格子里，退化成暴力遍历，所以格子尺寸最好设置成物体平均直径的 2 倍左右。

4.2 Narrow Phase：三种基本形状测试

Narrow Phase 是真正决定碰撞精确度的部分。我们的引擎有三个测试要写：圆-圆、圆-矩形（AABB）、矩形-矩形（AABB）。下面逐个拆解它们的数学原理。

圆与圆的检测最简单：如果两圆心距离小于半径之和，则相交。法线方向是连接圆心的方向，穿透深度是半径之和减去实际距离，接触点取两圆心连线的中点附近。注意不要直接开平方，先用length_sq和r_sum_sq比较，只在确认相交后才求平方根和法线。

def circle_circle_manifold(a, b): delta = b.position - a.position r_sum = a.shape_data.radius + b.shape_data.radius if delta.length_sq() > r_sum * r_sum: return None dist = delta.length() normal = delta * (1.0 / dist) if dist > 1e-8 else Vec2(1, 0) penetration = r_sum - dist contact_point = a.position + normal * (a.shape_data.radius - penetration * 0.5) return normal, penetration, contact_point

圆与矩形的检测思路是把圆心“夹”到矩形内部，得到一个最近点，然后看圆心到最近点的距离是否小于半径。这里的clamp操作是关键，它把圆心映射到矩形边界上最近的位置。如果圆心在矩形内部，最近点就等于圆心本身，这就需要一个退化分支：找出圆心离四条边最近的一面，以那面的法线作为碰撞法线。

def clamp_point_to_aabb(center, half_w, half_h): return Vec2( max(center.x - half_w, min(center.x, center.x)) if False else clamp(center.x, -half_w, half_w), clamp(center.y, -half_h, half_h) )

为了避免把代码写拧巴，我一般把 AABB 的位置字段视为“中心点”，矩形的表示是center + half_width + half_height，计算时一律先算绝对边界再与圆心比较。处理矩形内部的情况时，我会显式写一个“最近面”函数，把圆心相对于矩形中心的偏移分别投影到四条边的法线方向，取深度最小且为正的那个作为法线。这样代码更直白，也方便后续扩展旋转矩形。

矩形与矩形的 AABB 测试是整个引擎里最简单也最有价值的检测。两个 AABB 相交的条件是它们在两个轴向上的重叠都大于 0。重叠较少的轴，就是碰撞法线所在的轴，法线方向指向重叠为正的方向。

def aabb_aabb_manifold(a, b): dx = abs(b.position.x - a.position.x) dy = abs(b.position.y - a.position.y) ox = a.shape_data.half_w + b.shape_data.half_w - dx oy = a.shape_data.half_h + b.shape_data.half_h - dy if ox <= 0 or oy <= 0: return None if ox < oy: normal = Vec2(1, 0) if b.position.x > a.position.x else Vec2(-1, 0) penetration = ox contact_point = Vec2( (a.position.x + b.position.x) / 2.0, (a.position.y + b.position.y) / 2.0 ) else: normal = Vec2(0, 1) if b.position.y > a.position.y else Vec2(0, -1) penetration = oy contact_point = Vec2( (a.position.x + b.position.x) / 2.0, (a.position.y + b.position.y) / 2.0 ) return normal, penetration, contact_point

4.3 扩展到旋转矩形：SAT 分离轴定理的前瞻

如果哪一天你要支持矩形旋转，Narrow Phase 里的矩形-矩形检测就必须升级为分离轴定理（SAT）。SAT 的核心思想是：两个凸多边形不相交，当且仅当存在一条轴，使它们在这条轴上的投影不重叠；可能的分离轴恰好是每个多边形边的法线方向。因此你只需要取两个矩形总共四条边的法线，逐一做投影测试，就能判断是否相交，还能顺便求出最小穿透深度和法线。

在 AABB 版本里，检测轴正好就是 x 轴和 y 轴，这就是为什么前面那段代码看起来这么简洁。将来升级 OBB 时只需把“绝对边界比较”换成“把顶点坐标旋转”再投影。核心逻辑完全复用。我建议你在备注里写下这一点，方便以后翻代码时知道当初为什么要把 AABB 的先写出来。

4.4 接触流形：把几何信息打包传递给求解器

每次 Narrow Phase 得到的信息，我习惯用一个小结构体封装，就叫Contact：

class Contact: def __init__(self, body_a, body_b, normal, penetration, contact_point): self.a = body_a self.b = body_b self.normal = normal self.penetration = penetration self.point = contact_point

这里的关键点是：normal的方向约定必须统一。我统一约定为从 A 指向 B，这样冲量求解时只需要一个符号约定就能让两个物体分开。穿透深度是正值，接触点用来计算切向摩擦和冲量产生转动的力臂。把所有这几项打包后，后续的碰撞响应模块就可以完全不关心底层是圆还是矩形了——多态的效果。

5. 碰撞响应：冲量公式、恢复系数与摩擦力

碰撞检测告诉我们“物体撞上了”，碰撞响应则回答“接下来速度怎么变”。这里有一个基本物理概念上的分岔点：为什么用“冲量”而不是“力”。因为碰撞发生的时间极短，远小于一个模拟步长，力在这个过程中不是一个有限可测的量，真正起作用的是力对时间的积分，也就是冲量。它直接改变的是速度，而不是加速度。

5.1 法向冲量的公式推导

先假设两个物体只在接触点的法线方向上发生相互作用。定义法向相对速度vn = dot(vB - vA, normal)。如果vn > 0，说明两个物体正在分离，不需要响应；如果vn < 0，说明它们在相互接近，就必须施加冲量把速度改到分离方向。

设冲量大小为j，施加在 B 上的方向是normal，施加在 A 上的方向是-normal。速度变化量为：

ΔvA = -j * invMassA * normal ΔvB = j * invMassB * normal

碰撞后的法向相对速度应该满足恢复系数定义：

vn_after = -e * vn_before

把速度变化代入，展开后解出j，得到这个经典公式：

j = -(1 + e) * vn / (invMassA + invMassB)

注意分母是两者质量倒数之和，这正是为什么第一章节坚持用invMass的原因——公式里根本不需要单独写质量。如果某个物体是静态的，它的invMass = 0，分母自然只剩另一个物体的倒数，效果上等价于撞了一面无限重的墙。

落地代码如下。我把它写成一个独立的resolve_contact函数，方便后续进行迭代求解时反复调用。

def resolve_contact(contact, e): a, b = contact.a, contact.b n = contact.normal v_rel = b.velocity - a.velocity vn = v_rel.dot(n) if vn > 0.0: return # 已经分离 inv_mass_sum = a.inv_mass + b.inv_mass if inv_mass_sum < 1e-12: return # 两个物体都是静态 j = -(1.0 + e) * vn / inv_mass_sum impulse = n * j a.velocity -= impulse * a.inv_mass b.velocity += impulse * b.inv_mass

e是恢复系数，范围通常在 0 到 1 之间。e = 0代表完全非弹性，碰撞后法向相对速度为零，两个物体“粘”在一起；e = 1代表完全弹性，碰撞前后法向相对速度大小不变。真实世界里的物体都夹在中间，比如混凝土约等于 0.0，橡胶球约 0.8。如果你做的是休闲游戏，恢复系数给得低一些反而更容易稳定：高恢复系数会让能量积聚，堆叠场景很快会抖起来。

5.2 切向摩擦力：库仑模型的简化实现

法向冲量解决的是“弹开”，切向冲量解决的是“滑动”。库仑摩擦模型里，摩擦力的大小存在一个上限，与法向力成正比，这个比例系数就是摩擦系数 μ。摩擦系数分为静摩擦和动摩擦，简易引擎里通常只用一个“动态摩擦系数”近似。

计算切向相对速度的方法是：把相对速度向量减去其在法线上的投影分量，剩下的部分就是切向速度：

t = v_rel - n * (v_rel.dot(n)) if t.length_sq() > 1e-12: t = t.normalized() else: return # 没有切向运动 vt = v_rel.dot(t)

切向冲量的目标是把切向相对速度尽可能降低，但上限不能超过μ * |j|，这里的j是前面求出的法向冲量绝对值。这样做的好处是摩擦不会反向推动物体产生不自然的运动。实现代码如下：

def apply_friction(contact, mu, j_normal): a, b = contact.a, contact.b n = contact.normal v_rel = b.velocity - a.velocity t = v_rel - n * v_rel.dot(n) if t.length_sq() < 1e-12: return t = t.normalized() vt = v_rel.dot(t) inv_mass_sum = a.inv_mass + b.inv_mass jt = -vt / inv_mass_sum max_friction = mu * abs(j_normal) if jt > max_friction: jt = max_friction elif jt < -max_friction: jt = -max_friction impulse = t * jt a.velocity -= impulse * a.inv_mass b.velocity += impulse * b.inv_mass

mu取值的物理参考：冰面约 0.05，橡胶与路面约 1.0。在 Demo 里我会用 0.3 左右，既能看到滑行，也不会滑得收不住。

5.3 把法向冲量与摩擦合到一起的完整响应

完整响应函数就是先算法向冲量，再根据这个冲量计算摩擦。要注意顺序：法向冲量改变了速度，切向摩擦力基于的就是改变后的相对速度，因此两者在一个接触点上是连续计算，不是并行。把resolve_contact扩展成包含摩擦的版本，后面迭代求解的每一轮都调用它：

def solve_contact(contact, e, mu): a, b = contact.a, contact.b n = contact.normal v_rel = b.velocity - a.velocity vn = v_rel.dot(n) if vn > 0.0: return 0.0 inv_mass_sum = a.inv_mass + b.inv_mass j = -(1.0 + e) * vn / inv_mass_sum if inv_mass_sum > 1e-12 else 0.0 impulse = n * j a.velocity -= impulse * a.inv_mass b.velocity += impulse * b.inv_mass # 摩擦 v_rel2 = b.velocity - a.velocity t = v_rel2 - n * v_rel2.dot(n) if t.length_sq() > 1e-12: t = t.normalized() vt = v_rel2.dot(t) jt = -vt / inv_mass_sum if inv_mass_sum > 1e-12 else 0.0 max_f = mu * abs(j) if jt > max_f: jt = max_f elif jt < -max_f: jt = -max_f friction_impulse = t * jt a.velocity -= friction_impulse * a.inv_mass b.velocity += friction_impulse * b.inv_mass return j

到目前为止，所有的冲量计算都还停留在“两个物体之间的一次性碰撞”阶段。一旦物体数量变多，你会发现一个现象：底下的方块被压得不断下沉，上面的方块却还在乱弹。这就要进入下一章了。

6. 稳定堆叠：迭代求解与穿透修正

我敢说每个写过简易物理引擎的人，第一次把一堆方块叠起来时都经历过“快乐三秒然后全面崩塌”。原因是单个接触的冲量求解只顾局部，没有协调全局。打个比方：一个人站在一排人肩膀上，最底下的人受到的力并不是最顶上那个人直接压下来的力，而是整条压力链传导下来的总和，只求解一次接触，根本传递不到最底层。

6.1 为什么单次求解会让堆叠崩溃

假设三个方块 A、B、C 从上到下叠在一起，C 在最底，A 在最上。第一轮碰撞检测生成了两个接触对：A-B 对和 B-C 对。如果只对每个接触对求解一次，A 撞到 B 时 B 获得向下的速度增量，但 B 和 C 之间的接触在检测时发现 B 的速度已经是穿透方向，冲量会把它推回去一部分。这一套流程虽然能让每个局部碰撞看起来合理，但整体上 A 的重力势能还没有“传导”到底部，B 和 C 之间的接触力就不足以撑住 A，于是 B 慢慢下压，C 慢慢被压穿，几秒后全塌。

解决办法就是迭代。在同一帧内，对所有接触对反复求解多次，每次求解都会把上一轮产生的速度变化反馈到相邻接触，相当于在压力链上多次传递信息。Box2D 这类引擎里，那个velocityIterations参数就是干这个用的，常见取值在 8 到 20 之间。

迭代求解非常朴素，在一个循环里反复调用solve_contact即可：

def solve_all_contacts(contacts, iterations=10, e=0.3, mu=0.3): for _ in range(iterations): for contact in contacts: solve_contact(contact, e, mu)

第一次迭代后问题已经减轻，第二次之后压入深度会明显收敛，五六次后肉眼基本分辨不出堆叠下沉。代价是计算量线性增加，但对我们几十个物体的场景来说完全可以忽略。这个“简单到有点无脑”的方法，恰恰是工业级引擎里约束求解器的雏形——只不过工业级引擎加入了矩阵求解和更复杂的松弛策略。

6.2 穿透修正：为什么不能只用冲量

即使做了迭代，冲量求解也只能修正速度，已经发生的位移穿透并不会被完美拉回。想象一下，一个球陷入了地面 0.1 像素，冲量会让它的速度变为向上，但下一秒它仍然会带着那个负的位移重新落入地面，表现为肉眼可见的“地面之下抖动”。

常见的修正手段叫投影法（Projection）：直接把物体的位置沿法线方向推开穿透距离。但不能推满，否则物体会像被弹簧弹走一样弹出老远。我用的策略是：

• 设定slop（允许的微小穿透），低于这个值不修正，避免物体在稳定堆叠时一直微调抖动。slop通常取0.005。
• 设定修正比例percent，我习惯用0.8。一次修正 80% 的超出slop的穿透量，留下来的微小穿透让自然碰撞去缓慢归位。

def positional_correction(contacts, percent=0.8, slop=0.005): for contact in contacts: a, b = contact.a, contact.b inv_mass_sum = a.inv_mass + b.inv_mass if inv_mass_sum < 1e-12: continue correction_amount = max(contact.penetration - slop, 0.0) / inv_mass_sum correction = contact.normal * (correction_amount * percent) a.position -= correction * a.inv_mass b.position += correction * b.inv_mass

注意修正的是位置而不是速度，这本质上是在“作弊”——把物理上不该发生的穿透直接推回去。这个方法虽然不物理，但在游戏物理中完全够用，而且远比真正物理的“约束力求解”简单稳定。

6.3 迭代次数与主动休眠

迭代次数多到一定程度后效果增长非常有限，反而可能引入新的抖动。经验上 8 到 12 次已经是休闲游戏的安全区间；如果堆叠超过五层，可以考虑扩容到 20 次。省不省这笔计算，取决于你的目标平台——如果是移动端，每一层循环都在烧 CPU，值得用心调参。

另一个节省计算的手段是让静止物体进入休眠。休眠思想很简单：如果一个物体速度低于某个阈值（比如 0.01），且持续了若干帧，就停止更新它；只有当另一个活跃物体撞上它或者外力触发时才唤醒。休眠能显著减少堆叠场景的计算量，因为静止堆放弃后不再需要反复求解接触。简易引擎里可以先不做，但思路要留作扩展点。

6.4 从“抖”到“稳”的调参顺序总结

遇到堆叠不稳定，不要瞎调恢复系数和摩擦系数。先按这个顺序检查：

1. 物理步长是否固定？不固定，先修步长。
2. 迭代次数是否大于等于 8？没有，加上。
3. 穿透修正是否做了？没有，按percent=0.8, slop=0.005加上。
4. 恢复系数是否太高？建议先全部设为 0.1 ~ 0.3。
5. 摩擦系数是否太低？先试 0.5 再看效果。

这套顺序我踩坑踩出来的，照着调一遍，大部分抖动问题能解决八成以上。

7. 一个不加任何图形库的可运行 Demo

现在把所有模块拼起来，我写一个不依赖任何第三方库的完整脚本。它模拟了这样一幕：几只球和几个方块从半空掉下来，落在一片静止底板上，最终堆叠成一堆。每帧输出模拟数据，你可以直接跑，也可以把它封装进自己的渲染循环里。

这个脚本刻意不画画面，因为画出图形会引入平台相关的依赖，而物理引擎的核心逻辑是纯数值的。你完全可以用 HTML Canvas、Pygame、SDL 或者其它任何方式对接输出坐标数据。

7.1 完整代码：物理核心与模拟场景

下面这段代码总计不到两百行，就是一个可运行的最小物理引擎。我把前面所有讲解过的模块都整合进来了，注释里标明了每一段对应文章里的哪一章。

import random import math # ---------- Vec2 ---------- class Vec2: def __init__(self, x=0.0, y=0.0): self.x = float(x) self.y = float(y) def __add__(self, o): return Vec2(self.x + o.x, self.y + o.y) def __sub__(self, o): return Vec2(self.x - o.x, self.y - o.y) def __mul__(self, s): return Vec2(self.x * s, self.y * s) __rmul__ = __mul__ def dot(self, o): return self.x * o.x + self.y * o.y def cross(self, o): return self.x * o.y - self.y * o.x def length_sq(self): return self.x * self.x + self.y * self.y def length(self): return math.sqrt(self.length_sq()) def normalized(self): l = self.length() if l < 1e-10: return Vec2(1, 0) inv = 1.0 / l return Vec2(self.x * inv, self.y * inv) def perp(self): return Vec2(-self.y, self.x) # ---------- Body ---------- class Body: def __init__(self, shape, position, shape_data, mass=1.0, restitution=0.3): self.shape = shape # "circle" / "aabb" self.shape_data = shape_data # 半径 或 Vec2(半宽, 半高) self.position = position self.velocity = Vec2() self.angle = 0.0 self.angular_vel = 0.0 self.force_acc = Vec2() self.inv_mass = 0.0 if mass <= 0.0 else 1.0 / mass self.inv_i = 0.0 # 本Demo不使用旋转 self.restitution = restitution def apply_force(self, f): self.force_acc += f def step(self, gravity, dt): accel = self.force_acc * self.inv_mass self.velocity += accel * dt self.velocity += gravity * dt self.position += self.velocity * dt self.force_acc = Vec2() # ---------- 碰撞检测 ---------- def circle_circle(body_a, body_b): delta = body_b.position - body_a.position r_sum = body_a.shape_data + body_b.shape_data if delta.length_sq() > r_sum * r_sum: return None dist = delta.length() n = delta * (1.0 / dist) if dist > 1e-6 else Vec2(1, 0) pen = r_sum - dist cp = body_a.position + n * (body_a.shape_data - pen * 0.5) return n, pen, cp def aabb_aabb(body_a, body_b): a_hw = body_a.shape_data.x a_hh = body_a.shape_data.y b_hw = body_b.shape_data.x b_hh = body_b.shape_data.y dx = abs(body_b.position.x - body_a.position.x) dy = abs(body_b.position.y - body_a.position.y) ox = a_hw + b_hw - dx oy = a_hh + b_hh - dy if ox <= 0.0 or oy <= 0.0: return None if ox < oy: n = Vec2(1, 0) if body_b.position.x > body_a.position.x else Vec2(-1, 0) pen = ox else: n = Vec2(0, 1) if body_b.position.y > body_a.position.y else Vec2(0, -1) pen = oy cp = Vec2( (body_a.position.x + body_b.position.x) / 2.0, (body_a.position.y + body_b.position.y) / 2.0 ) return n, pen, cp def circle_aabb(body_circle, body_rect): rect_center = body_rect.position half_w = body_rect.shape_data.x half_h = body_rect.shape_data.y # 将圆心限制到矩形内部 closest = Vec2( max(rect_center.x - half_w, min(body_circle.position.x, rect_center.x + half_w)), max(rect_center.y - half_h, min(body_circle.position.y, rect_center.y + half_h)) ) delta = body_circle.position - closest if delta.length_sq() > body_circle.shape_data * body_circle.shape_data: return None if delta.length_sq() < 1e-12: # 圆心在矩形内部，用最近的面推出 overlap_left = body_circle.position.x - (rect_center.x - half_w) overlap_right = (rect_center.x + half_w) - body_circle.position.x overlap_bottom = body_circle.position.y - (rect_center.y - half_h) overlap_top = (rect_center.y + half_h) - body_circle.position.y min_ov = min(overlap_left, overlap_right, overlap_bottom, overlap_top) if min_ov == overlap_left: n = Vec2(-1, 0) elif min_ov == overlap_right: n = Vec2(1, 0) elif min_ov == overlap_bottom: n = Vec2(0, -1) else: n = Vec2(0, 1) pen = min_ov cp = body_circle.position + n * pen return n, pen, cp dist = delta.length() n = delta * (1.0 / dist) pen = body_circle.shape_data - dist cp = closest + n * (body_circle.shape_data - pen * 0.5) return n, pen, cp def detect_contacts(bodies): contacts = [] n = len(bodies) for i in range(n): for j in range(i + 1, n): a, b = bodies[i], bodies[j] if a.inv_mass == 0.0 and b.inv_mass == 0.0: continue if a.shape == "circle" and b.shape == "circle": result = circle_circle(a, b) elif a.shape == "aabb" and b.shape == "aabb": result = aabb_aabb(a, b) else: if a.shape == "circle": result = circle_aabb(a, b) else: # 交换顺序，让 circle_aabb 始终接收 (圆, 矩形)，法线取反 result = circle_aabb(b, a) if result: n_vec, pen, cp = result result = (n_vec * -1.0, pen, cp) if result: normal, penetration, cp = result # 注意：检测阶段不改变任何速度 contacts.append((a, b, normal, penetration, cp)) return contacts # ---------- 碰撞响应 ---------- def solve_contact(a, b, normal, penetration, cp, e, mu): v_rel = b.velocity - a.velocity vn = v_rel.dot(normal) if vn > 0.0: return 0.0 inv_sum = a.inv_mass + b.inv_mass if inv_sum < 1e-12: return 0.0 j = -(1.0 + e) * vn / inv_sum impulse = normal * j a.velocity -= impulse * a.inv_mass b.velocity += impulse * b.inv_mass # 摩擦 v_rel2 = b.velocity - a.velocity tangent = v_rel2 - normal * v_rel2.dot(normal) if tangent.length_sq() > 1e-10: tangent = tangent.normalized() vt = v_rel2.dot(tangent) jt = -vt / inv_sum max_friction = mu * abs(j) if jt > max_friction: jt = max_friction elif jt < -max_friction: jt = -max_friction f_impulse = tangent * jt a.velocity -= f_impulse * a.inv_mass b.velocity += f_impulse * b.inv_mass return j def solve_contacts(contacts, e, mu, iterations=10): for _ in range(iterations): for contact in contacts: solve_contact(contact[0], contact[1], contact[2], contact[3], contact[4], e, mu) def positional_correction(contacts, percent=0.8, slop=0.005): for contact in contacts: a, b = contact[0], contact[1] normal = contact[2] pen = contact[3] inv_sum = a.inv_mass + b.inv_mass if inv_sum < 1e-12: continue amount = max(pen - slop, 0.0) / inv_sum * percent corr = normal * amount a.position -= corr * a.inv_mass b.position += corr * b.inv_mass # ---------- 主模拟 ---------- def physics_step(bodies, gravity, dt, e=0.3, mu=0.3, iterations=10): # 1. 积分力与速度 for body in bodies: body.step(gravity, dt) # 2. 碰撞检测 contacts = detect_contacts(bodies) # 3. 迭代求解冲量与摩擦 solve_contacts(contacts, e, mu, iterations) # 4. 位置修正 positional_correction(contacts) def build_scene(): bodies = [] # 静态底板：一个巨大的AABB块 ground = Body("aabb", Vec2(0, -3.0), Vec2(10.0, 0.5), mass=0.0) bodies.append(ground) # 动态方块 for i in range(4): bodies.append(Body("aabb", Vec2(-1.0 + i * 0.7, 5.0 + i * 0.2), Vec2(0.35, 0.35), mass=1.0, restitution=0.1)) # 动态圆球 for i in range(4): bodies.append(Body("circle", Vec2(2.0, 4.0 + i * 0.5), 0.3, mass=1.0, restitution=0.3)) return bodies if __name__ == "__main__": random.seed(0) bodies = build_scene() gravity = Vec2(0.0, -9.81) dt = 1.0 / 60.0 for frame in range(180): physics_step(bodies, gravity, dt, e=0.2, mu=0.4, iterations=10) # 每30帧打印一次所有物体的位置 if frame % 30 == 0: info = [] for b in bodies: if b.inv_mass > 0.0: info.append(f"({b.position.x:.2f}, {b.position.y:.2f})") print(f"frame={frame:3d} objects: {' '.join(info)}")

把这段代码保存成.py文件直接运行，你会看到 180 帧（模拟 3 秒）内所有动态物体从半空掉落、碰撞、最终趋于静止的过程。如果某一次运行出现了堆叠下沉，检查是不是把迭代次数调得太低了；正常跑下来最后输出的一行，底部方块的位置应该稳定在y=2.15上下（静态地板中心在 -3.0，半高 0.5，顶面在 -2.5，方块半高 0.35，所以方块底面在 -2.5、中心应在 2.15）。这个数值就是静置状态的标尺。

7.2 怎样把数据变成可见画面

控制台输出只能验证引擎数值，但物理引擎的乐趣要看见才完整。如果你已经有一个渲染环境，把bodies里每个Body的位置和形状画出来即可。比如用 HTML Canvas，画一个圆就调用一次arc(x, y, r)，画一个矩形就调用rect(x - half_w, y - half_h, half_w * 2, half_h * 2)。没有任何其它魔法，物理核心永远只输出刚体状态数据。

如果你偏好 Python，那就在环境里安装 Pygame，写一个几十行的绘制循环；或者干脆把数据导出成 CSV，用任何电子表格软件画散点图。核心逻辑和画面完全解耦，这是这个设计最舒服的地方。

8. 再往前一步：从简易引擎到可商用引擎的路

写到这里，这套引擎已经能模拟基础刚体世界了，但它离一个“能扛事的引擎”还有相当距离。我最后一次打开自己的这个项目时，在代码旁边列了一张“还缺什么”的清单，应该对你有同样的参考价值。

8.1 下一步最值得做的四个扩展

扩展一：旋转刚体。给刚体加入角速度、角加速度，并把接触点偏移量写进冲量公式。核心变化是冲量分母多出两项：(rA × n)² * invIA + (rB × n)² * invIB，分别代表接触点冲量对两个物体转动的影响。扩展后，你的引擎就支持了最基础的刚体动力学，方块才能真正“翻滚”。

扩展二：SAT 支持任意凸多边形。轴对齐矩形版跑通后，换成 SAT 前文已经讲过原理，你可以慢慢把所有凸多边形统一到一个数据结构，碰撞检测会更通用，也能支持三角形。

扩展三：关节约束。距离关节、旋转关节、滑动关节，本质上是把“接触约束”泛化成“更一般的约束”。理解冲量法后，关节只是往约束矩阵里再塞几行约束方程。这类扩展是自制引擎最过瘾的部分，因为你会突然拥有“设计机关”的能力。

扩展四：Broad Phase 换成空间哈希。当物体数量超过几十个时，暴力遍历就开始拖后腿了。空间哈希每个格子里的物体维护一个链表，每帧只检测当前格子及邻居格子里的候选对。它的实现难度不大，但性能提升立竿见影。

8.2 用个人经验做个收尾，而不是替你做决定

说句心里话，现在常见的物理引擎都非常成熟，如果你的目标是做一款完整商业游戏，直接用成熟引擎是最快的路径。但我依然强烈建议每个想做物理相关功能的人亲手实现一遍简易引擎，哪怕只是跑通这个 Demo。原因不是“造轮子锻炼能力”这种空话，而是你在调堆叠参数、看它塌掉、再把它修好的过程中，会建立起对物理引擎内部工作的直觉。这个直觉非常值钱，它会在你以后用任何现成引擎时帮你快速定位问题，绕开那些让新人崩溃的坑。

如果你问我下一步怎么走，我会说：先把 Demo 加个简单的渲染层，让它动起来，感受一下自己亲手写出的“世界规则”。然后再给一个方块加上旋转，你要过第一个坎是冲量分母那两项怎么加——加对后，你会觉得整个世界都完整了。那一步的感受，跟看到小球第一次不穿透地面时一样爽。