4683 字

23 分钟

计算机图形学笔记(五)：动画与物理仿真

2025-07-24

计算机图形学

C++

/

编程

/

计算机图形学

/

OpenGL

NOTE
前四部分把”怎么把一帧画出来”讲完了。这部分换一个维度：怎么让一帧变成一段——动画与物理仿真。关键词是时间、插值、微分方程、约束。本章把骨骼蒙皮、弹簧质点、数值积分、约束求解的底层全串起来，为 Part 6 里的 PBD 布料、SPH 流体、FEM 做好数学准备。

目录#

动画基础 — 插值理论 / SLERP / Catmull-Rom 样条
关键帧与骨骼动画 — 关键帧曲线 / 骨骼层次 / LBS / DQS
牛顿力学回顾 — 位置-速度-加速度 / 牛顿第二定律 / 机械能
弹簧质点系统 — 胡克定律 / 阻尼 / 拓扑（绳/布）
数值积分 — 显式欧拉 / Verlet / RK4 / 隐式欧拉 / 辛 / 稳定性
约束求解 — 位置投影 / PBD / SHAKE
刚体与流体简介（Part 6 前导）

二十四、动画基础#

24.1 插值理论#

24.1.1 一维线性插值#

给定 $(t_0, f_0), (t_1, f_1)$ ，线性插值：

f(t) = (1-u) f_0 + u f_1, \qquad u = \frac{t - t_0}{t_1 - t_0} \in [0, 1]

性质：

端点插值: $f(t_0) = f_0, f(t_1) = f_1$
线性性:关于两端点值的加权平均
分量独立:对向量/颜色/多维值按分量各自插值

1
template <typename T>
2
T lerp(const T& a, const T& b, float u) {
3
    return a + u * (b - a);  // 一行版本，对浮点更稳定
4
}

TIP
写 a * (1-u) + b * u 在数值上不如 a + u * (b - a) 稳定，而且多一次乘法。

24.1.2 为什么线性插值对旋转不够#

用欧拉角线性插值会有万向节锁，用旋转矩阵逐元素插值会失去正交性。四元数是正确答案，但还要解决”单位球面上的插值”问题。

24.2 SLERP：球面线性插值#

给定单位四元数 $\mathbf{q}_0, \mathbf{q}_1$ ，夹角 $\theta = \arccos(\mathbf{q}_0 \cdot \mathbf{q}_1)$ :

\text{slerp}(\mathbf{q}_0, \mathbf{q}_1, u) = \frac{\sin((1-u)\theta)}{\sin\theta} \mathbf{q}_0 + \frac{\sin(u\theta)}{\sin\theta} \mathbf{q}_1

性质：

始终在单位球面上 → 结果仍是合法旋转
等角速度 → 动画看起来匀速旋转
最短大圆弧 → 走短路

📌 SLERP 的推导#

设 $\mathbf{q}(u) = a(u)\mathbf{q}_0 + b(u)\mathbf{q}_1$ 。三个条件：

$\mathbf{q}(0) = \mathbf{q}_0 \Rightarrow a(0) = 1, b(0) = 0$
$\mathbf{q}(1) = \mathbf{q}_1 \Rightarrow a(1) = 0, b(1) = 1$
$\|\mathbf{q}(u)\| = 1$ （单位性）

展开 $\|\mathbf{q}(u)\|^2 = a^2 + b^2 + 2ab\cos\theta = 1$ 。代入解析候选 $a = \sin((1-u)\theta)/\sin\theta, b = \sin(u\theta)/\sin\theta$ 验证即可。

1
Eigen::Quaternionf slerp(Eigen::Quaternionf q0, Eigen::Quaternionf q1, float u) {
2
    float cos_t = q0.dot(q1);
3
    // 选最短路径：q 与 -q 表示同一旋转
4
    if (cos_t < 0) { q1.coeffs() = -q1.coeffs(); cos_t = -cos_t; }
5
    // 极近时退化为 nlerp，避免 sin(theta) 趋零导致数值爆炸
6
    if (cos_t > 0.9995f) {
7
        Eigen::Quaternionf r(q0.coeffs() + u * (q1.coeffs() - q0.coeffs()));
8
        r.normalize();
9
        return r;
10
    }
11
    float theta = std::acos(cos_t);
12
    float s = std::sin(theta);
13
    float w0 = std::sin((1 - u) * theta) / s;
14
    float w1 = std::sin(u * theta) / s;
15
    return Eigen::Quaternionf(w0 * q0.coeffs() + w1 * q1.coeffs());
16
}

WARNING
引擎里骨骼旋转插值几乎全用 nlerp + 归一化 而非严格 SLERP——误差肉眼不可见但快 3 倍。SLERP 只在精度关键（长时间轨道、科学计算）时才用。

24.3 Catmull-Rom 样条（平滑曲线穿过控制点）#

Part 3 的贝塞尔/B 样条不一定经过控制点。动画里关键帧就是约束点，要求曲线必须穿过。Catmull-Rom 是经典选择：

\mathbf{C}(u) = \tfrac{1}{2}\begin{pmatrix} 1 & u & u^2 & u^3 \end{pmatrix} \begin{pmatrix} 0 & 2 & 0 & 0 \\ -1 & 0 & 1 & 0 \\ 2 & -5 & 4 & -1 \\ -1 & 3 & -3 & 1 \end{pmatrix} \begin{pmatrix} \mathbf{P}_0 \\ \mathbf{P}_1 \\ \mathbf{P}_2 \\ \mathbf{P}_3 \end{pmatrix}

曲线经过 $\mathbf{P}_1, \mathbf{P}_2$ ，切向量由 $(\mathbf{P}_2 - \mathbf{P}_0)/2$ 和 $(\mathbf{P}_3 - \mathbf{P}_1)/2$ 估计（中点差分）。 $C^1$ 连续，适合相机路径、UI 补间动画。

1
Eigen::Vector3f catmull_rom(const Eigen::Vector3f& p0, const Eigen::Vector3f& p1,
2
                            const Eigen::Vector3f& p2, const Eigen::Vector3f& p3, float u) {
3
    float u2 = u * u, u3 = u2 * u;
4
    return 0.5f * (
5
        (2.0f * p1) +
6
        (-p0 + p2) * u +
7
        (2.0f*p0 - 5.0f*p1 + 4.0f*p2 - p3) * u2 +
8
        (-p0 + 3.0f*p1 - 3.0f*p2 + p3) * u3
9
    );
10
}

二十五、关键帧与骨骼动画#

25.1 关键帧曲线#

关键帧 = 时间-值对（ $t_i, v_i$ ），中间时刻通过插值求值。工业引擎给每个关键帧配切线（线性/样条/阶跃/贝塞尔控制柄）来精细控制过渡。

1
template <typename T>
2
struct Keyframe {
3
    float t;
4
    T value;
5
    enum Interp { LINEAR, CUBIC, STEP } interp = LINEAR;
6
    T in_tangent, out_tangent;  // Hermite 插值用
7
};
8

9
template <typename T>
10
class AnimCurve {
11
    std::vector<Keyframe<T>> keys;   // 按 t 升序
12
public:
13
    T evaluate(float t) const {
14
        if (keys.empty()) return T{};
15
        if (t <= keys.front().t) return keys.front().value;
16
        if (t >= keys.back().t)  return keys.back().value;
17
        auto it = std::lower_bound(keys.begin(), keys.end(), t,
18
            [](const Keyframe<T>& k, float t){ return k.t < t; });
19
        const auto& b = *it;
20
        const auto& a = *(it - 1);
21
        float u = (t - a.t) / (b.t - a.t);
22
        switch (a.interp) {
23
            case Keyframe<T>::STEP:   return a.value;
24
            case Keyframe<T>::LINEAR: return lerp(a.value, b.value, u);
25
            case Keyframe<T>::CUBIC:  return hermite(a.value, a.out_tangent, b.value, b.in_tangent, u);
26
        }
27
        return a.value;
28
    }
29
};

25.2 骨骼层次#

骨骼是有根的树。每根骨头用两个变换刻画：

绑定姿势 bind pose 下的世界变换（静态）
当前姿势 current pose:每帧根据关键帧曲线+IK 算出的世界变换（动态）

顶点绑定时记录的是”相对于 bind pose 的位置”。渲染时先把顶点用 逆 bind 矩阵 拉回骨头局部空间，再用当前 world 矩阵推到新世界位置。

1
struct Bone {
2
    int parent = -1;
3
    Eigen::Affine3f local;           // 相对父骨的变换（关键帧驱动）
4
    Eigen::Affine3f world;           // 当前世界变换（每帧更新）
5
    Eigen::Affine3f inverse_bind;    // T-pose 下 world 的逆（bind 时算好）
6
};
7

8
void update_world_transforms(std::vector<Bone>& bones) {
9
    for (size_t i = 0; i < bones.size(); ++i) {
10
        const auto& b = bones[i];
11
        bones[i].world = (b.parent < 0)
12
            ? b.local
13
            : bones[b.parent].world * b.local;
14
    }
15
}

NOTE
父子顺序很关键：遍历时必须保证父节点已更新。工程上骨骼数组通常按拓扑序排列，正向一次扫过即可。

25.3 线性混合蒙皮（LBS, Linear Blend Skinning）#

每个顶点绑定若干骨头 + 权重（通常 ≤ 4）：

\mathbf{v}_{\text{skinned}} = \sum_{k} w_k \, (\mathbf{M}_k \cdot \mathbf{M}_k^{\text{bind}, -1}) \, \mathbf{v}_{\text{bind}}, \qquad \sum_k w_k = 1

1
Eigen::Vector3f skin_lbs(const Eigen::Vector3f& v_bind,
2
                         const std::array<int, 4>& bone_ids,
3
                         const std::array<float, 4>& weights,
4
                         const std::vector<Bone>& bones) {
5
    Eigen::Vector3f v_out = Eigen::Vector3f::Zero();
6
    for (int k = 0; k < 4; ++k) {
7
        if (weights[k] <= 0) continue;
8
        Eigen::Affine3f M = bones[bone_ids[k]].world * bones[bone_ids[k]].inverse_bind;
9
        v_out += weights[k] * (M * v_bind);
10
    }
11
    return v_out;
12
}

优点上只需加权矩阵 + 向量乘，实时渲染必备。

缺点:糖果扭曲（candy wrapper）——绕肘关节旋转 180° 时体积急剧塌缩。

25.4 双四元数蒙皮（DQS, Dual Quaternion Skinning）#

用对偶四元数表达刚体运动（旋转 + 平移），加权平均后仍保持等距——不会塌。Kavan 2008 论文是经典来源。

对偶四元数: $\hat{\mathbf{q}} = \mathbf{q}_r + \epsilon \mathbf{q}_d$ ，其中 $\mathbf{q}_r$ 是旋转部分， $\mathbf{q}_d = \tfrac{1}{2} \mathbf{t} \mathbf{q}_r$ 编码了平移。

📌 DQS 的加权流程#

每根骨头的世界变换转为对偶四元数 $\hat{\mathbf{q}}_k$
半球对齐:若 $\mathbf{q}_r^{(0)} \cdot \mathbf{q}_r^{(k)} < 0$ 则取反（和 SLERP 同样的理由）
加权和： $\hat{\mathbf{q}} = \sum_k w_k \hat{\mathbf{q}}_k$
归一化： $\hat{\mathbf{q}} \leftarrow \hat{\mathbf{q}} / \|\mathbf{q}_r\|$
作用在顶点上 $\mathbf{v}' = \hat{\mathbf{q}} \cdot \mathbf{v}$ （记住单位对偶四元数作用顶点的公式）

1
struct DualQuat {
2
    Eigen::Quaternionf real, dual;
3
    DualQuat operator+(const DualQuat& o) const {
4
        return { {real.coeffs() + o.real.coeffs()}, {dual.coeffs() + o.dual.coeffs()} };
5
    }
6
    DualQuat operator*(float s) const {
7
        return { {s * real.coeffs()}, {s * dual.coeffs()} };
8
    }
9
    void normalize() {
10
        float n = real.norm();
11
        real.coeffs() /= n;
12
        dual.coeffs() /= n;
13
    }
14
    Eigen::Vector3f transform(const Eigen::Vector3f& v) const {
15
        Eigen::Vector3f t = 2.0f * (dual * real.conjugate()).vec();
16
        return real * v + t;
17
    }
18
};

WARNING
DQS 对非均匀缩放无能为力（对偶四元数只表达刚体运动）。如果骨头有缩放，要么回退 LBS，要么做 LBS + DQS 混合蒙皮（Optimized Centers of Rotation Skinning 等现代算法）。

二十六、牛顿力学回顾#

NOTE
这一章是给后面弹簧质点、刚体、流体用的共同基石。符号与术语统一：位置 $\mathbf{x}$ 、速度 $\mathbf{v}$ 、加速度 $\mathbf{a}$ 、力 $\mathbf{F}$ 。

26.1 运动学#

\mathbf{v}(t) = \frac{d\mathbf{x}}{dt}, \qquad \mathbf{a}(t) = \frac{d\mathbf{v}}{dt} = \frac{d^2\mathbf{x}}{dt^2}

匀加速运动的解析解（ $\mathbf{a} = \text{const}$ ）:

\mathbf{v}(t) = \mathbf{v}_0 + \mathbf{a}t, \qquad \mathbf{x}(t) = \mathbf{x}_0 + \mathbf{v}_0 t + \tfrac{1}{2}\mathbf{a}t^2

26.2 牛顿第二定律#

\mathbf{F} = m\mathbf{a} = m\frac{d^2\mathbf{x}}{dt^2}

力的叠加:合力 = 各力的矢量和， $\mathbf{a} = (\sum_i \mathbf{F}_i)/m$ 。

26.3 常见力#

力	形式	备注
重力	$\mathbf{F}_g = m\mathbf{g}$	$\mathbf{g}$ 方向向下，量级 9.8 m/s² 或按场景缩放
弹簧力（胡克）	$\mathbf{F}_s = -k(\ell - L_0) \hat{\mathbf{e}}$	$\ell$ 当前长度， $L_0$ 原长， $\hat{\mathbf{e}}$ 沿弹簧单位向量
阻尼力	$\mathbf{F}_d = -\gamma \mathbf{v}$	线性阻尼， $\gamma$ 越大越耗能
库伦摩擦	$\\|\mathbf{F}_f\\| \leq \mu \\|\mathbf{F}_N\\|$	静摩擦 ≤ 阈值；动摩擦沿相对速度反向

26.4 机械能#

E = K + U, \qquad K = \tfrac{1}{2} m \|\mathbf{v}\|^2, \qquad U = mgh + \tfrac{1}{2} k (\ell - L_0)^2 + \cdots

能量漂移是评估数值积分器好坏的核心指标——保守系统的理想积分器应让 $E$ 长期保持有界（辛积分器做得到，显式欧拉做不到，见 §28）。

二十七、弹簧质点系统#

27.1 基本对象#

1
struct Mass {
2
    Eigen::Vector3f x;              // 当前位置
3
    Eigen::Vector3f x_prev;         // 上一步位置（Verlet 用）
4
    Eigen::Vector3f v;              // 速度
5
    Eigen::Vector3f f = Eigen::Vector3f::Zero();   // 累积力
6
    float m = 1.0f;
7
    bool pinned = false;            // 固定不动
8
};
9

10
struct Spring {
11
    int i, j;                        // 两端质点索引
12
    float k;                         // 弹簧系数
13
    float L0;                        // 原长
14
    float kd = 0.0f;                 // 结构阻尼（可选）
15
};

27.2 胡克力（含阻尼）#

带阻尼的弹簧力，沿弹簧方向的相对速度被阻尼掉：

\mathbf{F}_{ij} = -k(\ell - L_0)\hat{\mathbf{e}} - k_d (\mathbf{v}_{ij} \cdot \hat{\mathbf{e}}) \hat{\mathbf{e}}, \qquad \hat{\mathbf{e}} = \frac{\mathbf{x}_j - \mathbf{x}_i}{\ell}

1
void apply_spring(Mass& a, Mass& b, const Spring& s) {
2
    Eigen::Vector3f d = b.x - a.x;
3
    float len = d.norm();
4
    if (len < 1e-8f) return;
5
    Eigen::Vector3f e = d / len;
6
    Eigen::Vector3f v_rel = b.v - a.v;
7
    Eigen::Vector3f F = s.k * (len - s.L0) * e
8
                      + s.kd * v_rel.dot(e) * e;
9
    a.f +=  F;       // 牛顿第三定律
10
    b.f += -F;
11
}

27.3 典型拓扑#

27.3.1 绳（1D 链）#

线性串接质点：springs.emplace_back(i, i+1, k, L0)，一端或两端 pin 住。

27.3.2 布料（2D 网格）#

弹簧种类	连接方式	物理作用
结构弹簧	水平 + 垂直邻居	抵抗拉伸
剪切弹簧	对角线	抵抗面内剪切形变
弯曲弹簧	跨一个节点的直线连接	抵抗面外弯折

1
void build_cloth_springs(int rows, int cols, float k_struct, float k_shear, float k_bend,
2
                         float dx, std::vector<Spring>& out) {
3
    auto id = [cols](int r, int c){ return r * cols + c; };
4
    // 结构
5
    for (int r = 0; r < rows; ++r)
6
        for (int c = 0; c < cols; ++c) {
7
            if (c + 1 < cols) out.push_back({id(r,c), id(r,c+1), k_struct, dx});
8
            if (r + 1 < rows) out.push_back({id(r,c), id(r+1,c), k_struct, dx});
9
        }
10
    // 剪切
11
    for (int r = 0; r + 1 < rows; ++r)
12
        for (int c = 0; c + 1 < cols; ++c) {
13
            float d = dx * std::sqrt(2.0f);
14
            out.push_back({id(r,c),   id(r+1,c+1), k_shear, d});
15
            out.push_back({id(r+1,c), id(r,c+1),   k_shear, d});
16
        }
17
    // 弯曲
18
    for (int r = 0; r < rows; ++r)
19
        for (int c = 0; c + 2 < cols; ++c)
20
            out.push_back({id(r,c), id(r,c+2), k_bend, 2*dx});
21
    for (int r = 0; r + 2 < rows; ++r)
22
        for (int c = 0; c < cols; ++c)
23
            out.push_back({id(r,c), id(r+2,c), k_bend, 2*dx});
24
}

27.3.3 杨氏模量与 k 的关系#

弹簧系数 $k$ 可以由材料参数推出：

k = \frac{E \cdot A}{L_0}

$E$ :杨氏模量（Pa）
$A$ :等效截面积
$L_0$ :原长

这让”弹簧系数”从凭感觉调变成”给材料类型”。

二十八、数值积分（核心）#

WARNING
整个物理仿真的工程难点几乎都压在数值积分器身上。选错会发散、能量漂移、抖动。先记住这三条原则：刚度越高，步长越小、显式方法能量正漂、隐式方法能量负漂、保守系统用辛积分器。

28.1 显式欧拉#

\mathbf{v}^{n+1} = \mathbf{v}^n + \Delta t \, \mathbf{a}^n, \qquad \mathbf{x}^{n+1} = \mathbf{x}^n + \Delta t \, \mathbf{v}^n

优点:最简单，一眼看懂。

缺点:条件稳定 —— 对弹簧系统要求 $\Delta t < 2/\omega, \omega = \sqrt{k/m}$ ，且总是把能量注入系统（正能量漂移），长时间仿真会爆炸。

28.2 辛欧拉（Symplectic Euler，半隐式）#

把位置更新改用新速度:

\mathbf{v}^{n+1} = \mathbf{v}^n + \Delta t \, \mathbf{a}^n, \qquad \mathbf{x}^{n+1} = \mathbf{x}^n + \Delta t \, \mathbf{v}^{n+1}

代码改动只有一行，但能量长期有界，游戏物理里的默认选择。

1
void symplectic_euler(Mass& p, float dt) {
2
    if (p.pinned) { p.f.setZero(); return; }
3
    p.v += dt * (p.f / p.m);
4
    p.x += dt * p.v;
5
    p.f.setZero();
6
}

28.3 Verlet 积分#

28.3.1 位置 Verlet#

消掉速度，只看位置：

\mathbf{x}^{n+1} = 2\mathbf{x}^n - \mathbf{x}^{n-1} + \Delta t^2 \, \mathbf{a}^n

速度作为差分估计 $\mathbf{v} \approx (\mathbf{x}^{n+1} - \mathbf{x}^{n-1})/(2\Delta t)$ 。约束投影型仿真（PBD、布料）首选。

1
void position_verlet(Mass& p, float dt) {
2
    if (p.pinned) { p.f.setZero(); return; }
3
    Eigen::Vector3f tmp = p.x;
4
    Eigen::Vector3f a = p.f / p.m;
5
    p.x = 2.0f * p.x - p.x_prev + a * dt * dt;
6
    p.x_prev = tmp;
7
    p.f.setZero();
8
}

28.3.2 速度 Verlet#

\mathbf{x}^{n+1} = \mathbf{x}^n + \Delta t \, \mathbf{v}^n + \tfrac{1}{2}\Delta t^2 \, \mathbf{a}^n

\mathbf{v}^{n+1} = \mathbf{v}^n + \tfrac{1}{2}\Delta t \, (\mathbf{a}^n + \mathbf{a}^{n+1})

二阶精度 + 辛，做分子动力学与天体仿真几乎必选。

28.4 Runge-Kutta 4（RK4）#

经典四阶显式法，精度远高于欧拉：

\begin{aligned} \mathbf{k}_1 &= f(\mathbf{x}^n,\, t^n), \\ \mathbf{k}_2 &= f(\mathbf{x}^n + \tfrac{\Delta t}{2}\mathbf{k}_1,\, t^n + \tfrac{\Delta t}{2}), \\ \mathbf{k}_3 &= f(\mathbf{x}^n + \tfrac{\Delta t}{2}\mathbf{k}_2,\, t^n + \tfrac{\Delta t}{2}), \\ \mathbf{k}_4 &= f(\mathbf{x}^n + \Delta t\,\mathbf{k}_3,\, t^n + \Delta t), \\ \mathbf{x}^{n+1} &= \mathbf{x}^n + \tfrac{\Delta t}{6}(\mathbf{k}_1 + 2\mathbf{k}_2 + 2\mathbf{k}_3 + \mathbf{k}_4). \end{aligned}

WARNING
RK4 不是辛的——对保守系统仍会缓慢漂移能量。高刚度系统用 RK4 仍可能不稳（刚性问题 → 用隐式法）。

28.5 隐式欧拉#

\mathbf{v}^{n+1} = \mathbf{v}^n + \Delta t \, \mathbf{a}(\mathbf{x}^{n+1}), \qquad \mathbf{x}^{n+1} = \mathbf{x}^n + \Delta t \, \mathbf{v}^{n+1}

右侧的 $\mathbf{a}(\mathbf{x}^{n+1})$ 依赖未知的 $\mathbf{x}^{n+1}$ ，要解一个线性/非线性方程组。Baraff-Witkin 1998 的布料仿真就是隐式积分的开山之作——线性化后解一个稀疏线性系统：

(\mathbf{M} - \Delta t^2 \mathbf{K} - \Delta t \mathbf{D}) \Delta \mathbf{v} = \Delta t (\mathbf{F} + \Delta t \mathbf{K} \mathbf{v})

其中 $\mathbf{K}$ 是力的雅可比（刚度矩阵）， $\mathbf{D}$ 是阻尼雅可比。

优点:无条件稳定。

代价:每步都要解线性系统；数值阻尼 使能量长期衰减（这正是”布料运动看起来黏黏的”的来源）。

28.6 稳定性与 CFL 条件#

28.6.1 弹簧-质点系统的显式稳定性#

线性振子 $\ddot{x} = -\omega^2 x$ ， $\omega = \sqrt{k/m}$ 。显式欧拉的放大因子 $|1 + i\omega\Delta t|$ 总 > 1 → 任意步长都不稳（严格说能量单调增）。辛欧拉稳定条件：

\Delta t < \frac{2}{\omega} = 2\sqrt{\frac{m}{k}}

28.6.2 CFL 条件（波方程）#

信息传播速度 $c$ ，空间离散 $\Delta x$ ，则：

\Delta t \leq \alpha \frac{\Delta x}{c}, \qquad \alpha \approx 0.5

对弹簧： $c = \sqrt{k / m_{\text{reduced}}}$ ， $\Delta x = L_0$ ，CFL 给出与稳定性条件一致的上限。实用建议：在稳定性阈值的 50% 以下选步长，留一半余量给非线性。

二十九、约束求解#

WARNING
很多”物理现象”本质是约束:绳不可拉伸、刚体两点距离不变、角色脚不穿地。直接用胡克建模就得要超大 $k$ ，逼到显式积分步长小得跑不动。把约束当约束本身求解（而不是当作软弹簧）是现代实时物理的主线。

29.1 位置投影#

距离约束 $C(\mathbf{x}_i, \mathbf{x}_j) = \|\mathbf{x}_j - \mathbf{x}_i\| - L_0 = 0$ 。任意时刻若偏离，做一次投影修正:

\Delta \mathbf{x}_i = + \frac{w_i}{w_i + w_j} C \hat{\mathbf{e}}, \qquad \Delta \mathbf{x}_j = -\frac{w_j}{w_i + w_j} C \hat{\mathbf{e}}

$w_i = 1/m_i$ （固定点 $w = 0$ ）
$\hat{\mathbf{e}}$ 为沿当前连线的单位向量
按反比质量分配修正,保证动量守恒

1
void project_distance(Mass& a, Mass& b, float L0) {
2
    float wa = a.pinned ? 0.0f : 1.0f / a.m;
3
    float wb = b.pinned ? 0.0f : 1.0f / b.m;
4
    float w_sum = wa + wb;
5
    if (w_sum < 1e-8f) return;
6
    Eigen::Vector3f d = b.x - a.x;
7
    float len = d.norm();
8
    if (len < 1e-8f) return;
9
    Eigen::Vector3f e = d / len;
10
    float C = len - L0;
11
    a.x += (wa / w_sum) * C * e;
12
    b.x -= (wb / w_sum) * C * e;
13
}

29.2 Position Based Dynamics (PBD)#

Müller 2006/2007 提出的现代范式。核心套路：

用显式/辛欧拉预测 $\tilde{\mathbf{x}}$
迭代投影每条约束到 $\tilde{\mathbf{x}}$ 上
从修正后的位置反推新速度 $\mathbf{v}^{n+1} = (\mathbf{x}^{n+1} - \mathbf{x}^n)/\Delta t$

1
void pbd_step(std::vector<Mass>& ps, const std::vector<Spring>& cs,
2
              const Eigen::Vector3f& g, float dt, int iters) {
3
    // 1) 预测
4
    std::vector<Eigen::Vector3f> x_pred(ps.size());
5
    for (size_t i = 0; i < ps.size(); ++i) {
6
        if (ps[i].pinned) { x_pred[i] = ps[i].x; continue; }
7
        ps[i].v += dt * g;
8
        x_pred[i] = ps[i].x + dt * ps[i].v;
9
    }
10
    // 2) 约束投影
11
    for (int it = 0; it < iters; ++it) {
12
        for (const auto& s : cs) {
13
            Eigen::Vector3f& xi = x_pred[s.i];
14
            Eigen::Vector3f& xj = x_pred[s.j];
15
            Eigen::Vector3f d = xj - xi;
16
            float len = d.norm();
17
            if (len < 1e-8f) continue;
18
            Eigen::Vector3f e = d / len;
19
            float C = len - s.L0;
20
            float wi = ps[s.i].pinned ? 0.0f : 1.0f / ps[s.i].m;
21
            float wj = ps[s.j].pinned ? 0.0f : 1.0f / ps[s.j].m;
22
            float w_sum = wi + wj;
23
            if (w_sum < 1e-8f) continue;
24
            xi += (wi / w_sum) * C * e;
25
            xj -= (wj / w_sum) * C * e;
26
        }
27
    }
28
    // 3) 更新速度与位置
29
    for (size_t i = 0; i < ps.size(); ++i) {
30
        if (ps[i].pinned) continue;
31
        ps[i].v = (x_pred[i] - ps[i].x) / dt;
32
        ps[i].x = x_pred[i];
33
    }
34
}

PBD 的性格:

几乎无条件稳定（都是位置修正，不会像力爆炸）
刚度通过”迭代次数”和”每次投影的比例因子”间接控制——不像胡克 $k$ 直接对应物理
XPBD (2016) 补上了拉格朗日乘子的严谨版本，让刚度变成物理参数

TIP
游戏里 90% 的布料、绳、毛发都是 PBD/XPBD 实现。它和 Part 6 的 FEM/SPH 并不互斥——PBD 可以当作 FEM 的快速近似。

29.3 SHAKE / RATTLE（分子动力学侧路线）#

SHAKE 是 Verlet 上的约束扩展，用拉格朗日乘子迭代直到约束误差小于阈值：

C(\mathbf{x}^{n+1}) = 0 \Rightarrow \mathbf{x}^{n+1} \leftarrow \mathbf{x}^{n+1} - \lambda \nabla C / m

比 PBD 更强调数学严格性（蒙卡尔洛分子动力学要求严格保约束），图形里用得少，但理解它就知道 PBD 是一种”用位置投影逼近拉格朗日乘子”的近亲。

三十、刚体与流体概述（Part 6 前导）#

30.1 刚体#

质点的推广：带旋转的状态 $(\mathbf{x}, \mathbf{v}, \mathbf{q}, \boldsymbol{\omega})$ 。核心方程：

m\dot{\mathbf{v}} = \mathbf{F}, \qquad \mathbf{I}(t)\dot{\boldsymbol{\omega}} + \boldsymbol{\omega} \times (\mathbf{I}(t)\boldsymbol{\omega}) = \boldsymbol{\tau}

$\mathbf{I}(t) = \mathbf{R}(t) \mathbf{I}_{\text{body}} \mathbf{R}(t)^T$ :世界坐标惯性张量
$\boldsymbol{\tau}$ :合力矩
四元数运动学： $\dot{\mathbf{q}} = \tfrac{1}{2} \boldsymbol{\omega} \mathbf{q}$ （把角速度当纯四元数左乘）

碰撞响应用冲量法 $\Delta \mathbf{v} = j \hat{\mathbf{n}} / m$ ， $j$ 由恢复系数 + 库伦摩擦锥求解。Part 6 会展开。

30.2 流体（欧拉 vs 拉格朗日）#

欧拉法（网格） 2003 “Stable Fluids”，在固定网格上解 Navier-Stokes，烟、火焰常用
拉格朗日法（粒子） (Smoothed Particle Hydrodynamics)、PBF (Position Based Fluids)，液体自由表面常用

两者的核心都是 Navier-Stokes：

\frac{\partial \mathbf{u}}{\partial t} + (\mathbf{u} \cdot \nabla)\mathbf{u} = -\frac{1}{\rho}\nabla p + \nu \nabla^2 \mathbf{u} + \mathbf{f}

完整推导、SPH 核函数、PCISPH/PBF 的求解流程放 Part 6。

30.3 Part 5 → Part 6 的进度条#

主题	本部分做到	Part 6 展开
弹簧质点	推导 + 实现 + 拓扑	XPBD、FEM、有限元连续体
积分器	欧拉 / Verlet / RK4 / 隐式 / 辛	指数积分器、IMEX 多尺度
约束	PBD 基本框架	XPBD 数学严格、收缩映射证明
刚体	牛-欧拉方程 + 惯性张量	冲量法、LCP 约束接触
流体	Navier-Stokes 公式	Stable Fluids / SPH / PBF

小结#

主题	工具	canonical 位置
标量/向量插值	lerp, Hermite, Catmull-Rom	§24
旋转插值	SLERP / nlerp	§24.2 — 骨骼动画的基础
骨骼与蒙皮	骨骼层次 / LBS / DQS	§25
经典力学	牛顿第二定律 / 机械能	§26 — 后面所有仿真共用的底
弹簧质点	胡克 + 阻尼 + 结构/剪切/弯曲	§27
数值积分	显式/辛/Verlet/RK4/隐式	§28 — 稳定性决定工程选择
约束投影	PBD / SHAKE / 距离约束	§29 — 实时布料与绳
刚体&流体	牛-欧拉方程 / Navier-Stokes	§30 概述 → Part 6 展开