XPBD公式推导：从隐式欧拉到可执行算法

这是一份从零开始的XPBD（Extended Position Based Dynamics）完整推导文档。我们将从最基础的物理原理出发，通过清晰的逻辑链条，最终得到可在计算机上实现的算法。

---

第一部分：基础概念

1.1 牛顿第二定律

概念：物体的加速度与所受合外力成正比，与质量成反比。

公式：

$$\mathbf{F}=m\mathbf{a}=m\ddot{\mathbf{p}}$$

例子：一个1kg的小球，受到10N的力，加速度为10 m/s²。

1.2 时间离散化

概念：计算机无法处理连续时间，需要将时间分成小段（时间步长 $\Delta t$）。

例子：$\Delta t=0.01$秒表示每秒更新100次。

1.3 欧拉积分

概念：用差分近似微分，从当前状态推算下一时刻状态。

前向欧拉（显式）：

$${\mathbf{v}}^{n+1}={\mathbf{v}}^n+\Delta t\cdot {\mathbf{M}}^{-1}\mathbf{F}({\mathbf{p}}^n)$$ $${\mathbf{p}}^{n+1}={\mathbf{p}}^n+\Delta t\cdot {\mathbf{v}}^{n+1}$$

问题：不稳定，大步长容易“爆炸”。

1.4 隐式欧拉

概念：使用未来时刻的力来计算当前更新，数值稳定。

公式：

$${\mathbf{v}}^{n+1}={\mathbf{v}}^n+\Delta t\cdot {\mathbf{M}}^{-1}\mathbf{F}({\mathbf{p}}^{n+1})$$ $${\mathbf{p}}^{n+1}={\mathbf{p}}^n+\Delta t\cdot {\mathbf{v}}^{n+1}$$

问题：${\mathbf{p}}^{n+1}$ 同时出现在两边，需要解方程。

例子：弹簧力 $F=-kx$，隐式欧拉会产生一个关于 $x^{n+1}$ 的方程：

$$x^{n+1}=x^n+\Delta t{v}^n-\frac{k}{m}\Delta t^2{x}^{n+1}$$

---

第二部分：从隐式欧拉到优化问题

2.1 消去速度

将第二式代入第一式：

$$\frac{{\mathbf{p}}^{n+1}-{\mathbf{p}}^n}{\Delta t}={\mathbf{v}}^n+\Delta t\cdot {\mathbf{M}}^{-1}\mathbf{F}({\mathbf{p}}^{n+1})$$

2.2 引入预测位置

定义预测位置：

$${\mathbf{p}}^{\ast }={\mathbf{p}}^n+\Delta t{\mathbf{v}}^n$$

代入得：

$$\boxed{{\mathbf{p}}^{n+1}-{\mathbf{p}}^{\ast }=\Delta t^2{\mathbf{M}}^{-1}\mathbf{F}({\mathbf{p}}^{n+1})}\text{\hspace{1em}(1)}$$

物理意义：位置修正量 = 时间步长平方 × 加速度。

2.3 保守力与势能

概念：保守力是可以用势能函数描述的力，满足 $\mathbf{F}=-\nabla U$。

例子：重力势能 $U=mgh$，重力 $\mathbf{F}=(0,-mg)$。

在XPBD中，约束力来自约束势能：

$$U_{\text{con}}=\frac{1}{2}{\alpha }^{-1}C(\mathbf{p}{)}^2$$

其中 $C(\mathbf{p})$ 是约束函数（如距离约束 $|{\mathbf{p}}_2-{\mathbf{p}}_1|-d$）。

2.4 转化为优化问题

关键定理：方程 $\mathbf{p}-{\mathbf{p}}^{\ast }=\Delta t^2{\mathbf{M}}^{-1}\mathbf{F}(\mathbf{p})$ 的解，等价于以下最小化问题的解：

$$\boxed{\underset{\mathbf{p}}{\min }\Phi (\mathbf{p})=\underset{\mathbf{p}}{\min }\left[\frac{1}{2\Delta t^2}\Vert \mathbf{p}-{\mathbf{p}}^{\ast }{\Vert }_{\mathbf{M}}^2+U(\mathbf{p})\right]}\text{\hspace{1em}(2)}$$

验证：对 $\Phi$ 求梯度：

$$\nabla \Phi (\mathbf{p})=\frac{1}{\Delta t^2}\mathbf{M}(\mathbf{p}-{\mathbf{p}}^{\ast })+\nabla U(\mathbf{p})=\frac{1}{\Delta t^2}\mathbf{M}(\mathbf{p}-{\mathbf{p}}^{\ast })-\mathbf{F}(\mathbf{p})$$

令梯度为零，即得方程(1)。

---

第三部分：Legendre变换

3.1 为什么需要Legendre变换？

问题：式(2)中的势能项 $\frac{1}{2}{\alpha }^{-1}C(\mathbf{p}{)}^2$ 是 $C(\mathbf{p})$ 的二次函数，而 $C(\mathbf{p})$ 本身是非线性的。直接处理很困难。

目标：引入一个辅助变量 $\lambda$，将二次项变成线性项。

3.2 Legendre变换的定义

概念：对于凸函数 $f(x)=\frac{1}{2}a{x}^2$，其Legendre变换为：

$$f^{\ast }(p)=\underset{p}{\max }[px-f(x)]$$

推导：

$$\frac{d}{dx}(px-\frac{1}{2}a{x}^2)=p-ax=0\Rightarrow x=p/a$$ $$f^{\ast }(p)=p\cdot (p/a)-\frac{1}{2}a(p/a{)}^2=\frac{p^2}{a}-\frac{p^2}{2a}=\frac{p^2}{2a}$$

恒等式：

$$\boxed{\frac{1}{2}a{x}^2=\underset{p}{\max }\left(px-\frac{1}{2a}{p}^2\right)}\text{\hspace{1em}(3)}$$

例子：$a=2,x=3$，左边 $=9$；右边对 $p$ 最大化，$p=ax=6$，$6\times 3-\frac{1}{4}\times 36=18-9=9$。

3.3 应用到约束势能

令 $a={\alpha }^{-1}$，$x=C(\mathbf{p})$，$p=\lambda$：

$$\boxed{\frac{1}{2}{\alpha }^{-1}C(\mathbf{p}{)}^2=\underset{\lambda }{\max }\left(\lambda C(\mathbf{p})-\frac{1}{2}\alpha {\lambda }^2\right)}\text{\hspace{1em}(4)}$$

3.4 重写优化问题

将(4)代入(2)：

$$\underset{\mathbf{p}}{\min }\underset{\mathbf{\lambda }}{\max }\Psi (\mathbf{p},\mathbf{\lambda })=\underset{\mathbf{p}}{\min }\underset{\mathbf{\lambda }}{\max }\left[\frac{1}{2\Delta t^2}\Vert \mathbf{p}-{\mathbf{p}}^{\ast }{\Vert }_{\mathbf{M}}^2+\sum _j\left({\lambda }_j{C}_j(\mathbf{p})-\frac{1}{2}{\alpha }_j{\lambda }_j^2\right)\right]\text{\hspace{1em}(5)}$$

这是一个鞍点问题：对 $\mathbf{p}$ 最小化，对 $\lambda$ 最大化。

---

第四部分：最优性条件

4.1 对位置求导

固定 $\lambda$，对 ${\mathbf{p}}_i$ 求导：

$$\frac{\partial \Psi }{\partial {\mathbf{p}}_i}=\frac{m_i}{\Delta t^2}({\mathbf{p}}_i-{\mathbf{p}}_i^{\ast })+\sum _j{\lambda }_j{\nabla }_{{\mathbf{p}}_i}{C}_j=0$$

得到位置更新公式：

$$\boxed{{\mathbf{p}}_i={\mathbf{p}}_i^{\ast }-\frac{\Delta t^2}{m_i}\sum _j{\lambda }_j{\nabla }_{{\mathbf{p}}_i}{C}_j}\text{\hspace{1em}(6)}$$

4.2 对 $\lambda$ 求导

固定 $\mathbf{p}$，对 ${\lambda }_j$ 求导：

$$\frac{\partial \Psi }{\partial {\lambda }_j}=C_j(\mathbf{p})-{\alpha }_j{\lambda }_j=0$$

得到约束条件：

$$\boxed{C_j(\mathbf{p})={\alpha }_j{\lambda }_j}\text{\hspace{1em}(7)}$$

---

第五部分：线性化与单约束求解

5.1 为什么需要线性化？

问题：式(6)和(7)是耦合的非线性方程组，不能直接求解。

策略：在当前位置 ${\mathbf{p}}^0$ 附近线性化约束，然后一次只处理一个约束。

5.2 泰勒展开

$$C_j(\mathbf{p})\approx C_j({\mathbf{p}}^0)+\sum _i{\nabla }_{{\mathbf{p}}_i}{C}_j\cdot ({\mathbf{p}}_i-{\mathbf{p}}_i^0)\text{\hspace{1em}(8)}$$

5.3 代入位置更新

对于单个约束 $j$，假设其他约束的贡献已被吸收进当前位置，我们有：

$${\mathbf{p}}_i-{\mathbf{p}}_i^0=-\frac{\Delta t^2}{m_i}{\lambda }_j^{\text{new}}{\nabla }_{{\mathbf{p}}_i}{C}_j$$

代入(8)并利用(7)：

$$C_j({\mathbf{p}}^0)-\Delta t^2{\lambda }_j^{\text{new}}\sum _i\frac{1}{m_i}\Vert {\nabla }_{{\mathbf{p}}_i}{C}_j{\Vert }^2={\alpha }_j{\lambda }_j^{\text{new}}$$

记 ${\text{denom}}_j={\sum }_i{w}_i\Vert {\nabla }_{{\mathbf{p}}_i}{C}_j{\Vert }^2$，其中 $w_i=1/m_i$。

5.4 解出 ${\lambda }_j^{\text{new}}$

$${\lambda }_j^{\text{new}}=\frac{C_j({\mathbf{p}}^0)}{{\alpha }_j+\Delta t^2\cdot {\text{denom}}_j}\text{\hspace{1em}(9)}$$

5.5 增量形式

设 ${\lambda }_j^{\text{new}}={\lambda }_j^{\text{old}}+\Delta {\lambda }_j$，代入(9)：

$$\Delta {\lambda }_j=\frac{C_j({\mathbf{p}}^0)-{\alpha }_j{\lambda }_j^{\text{old}}}{{\alpha }_j+\Delta t^2\cdot {\text{denom}}_j}\text{\hspace{1em}(10)}$$

物理意义：分子是当前约束误差减去已有约束力产生的形变。

---

第六部分：Gauss-Seidel迭代

6.1 什么是Gauss-Seidel迭代？

概念：求解方程组的一种方法，每次更新一个变量，立即使用最新值。

例子：解方程组

$$\{\begin{array}{l}x+2y=5\\ 3x+y=6\end{array}$$

Gauss-Seidel迭代：

1. 假设 $y=0$，得 $x=5$
2. 用 $x=5$，得 $3\times 5+y=6\Rightarrow y=-9$
3. 用 $y=-9$，得 $x-18=5\Rightarrow x=23$
4. 收敛到 $(1.4,1.8)$

6.2 为什么需要Gauss-Seidel？

问题：多个约束共享粒子，不能独立求解。

XPBD的策略：

• 一次只处理一个约束
• 立即更新受影响的粒子位置
• 下一个约束看到的是更新后的位置
• 多次迭代逼近全局解

6.3 完整算法

输入: 当前位置 p, 逆质量 w, 时间步长 Δt
输出: 更新后的位置 p
 
for iteration = 1 to solver_iterations:
    for each constraint j:
        # 计算当前约束值
        C = C_j(p)
        
        # 计算梯度平方和
        denom = Σ_i w_i * ||∇C_i||²
        
        # 更新 λ
        Δλ = (C - α_j * λ_j) / (α_j + Δt² * denom)
        λ_j = λ_j + Δλ
        
        # 更新位置
        for each particle i in constraint j:
            Δp = -w_i * Δt² * Δλ * ∇C_i
            p_i = p_i + Δp

---

第七部分：距离约束的具体推导

7.1 约束函数

$$C({\mathbf{p}}_1,{\mathbf{p}}_2)=|{\mathbf{p}}_2-{\mathbf{p}}_1|-d$$

7.2 梯度计算

设 $\mathbf{n}=\frac{{\mathbf{p}}_2-{\mathbf{p}}_1}{|{\mathbf{p}}_2-{\mathbf{p}}_1|}$，则：

$${\nabla }_{{\mathbf{p}}_1}C=-\mathbf{n},{\nabla }_{{\mathbf{p}}_2}C=\mathbf{n}$$

7.3 分母计算

$$\text{denom}=w_1\Vert -\mathbf{n}{\Vert }^2+w_2\Vert \mathbf{n}{\Vert }^2=w_1+w_2$$

7.4 修正公式

$$\Delta \lambda =\frac{|{\mathbf{p}}_2-{\mathbf{p}}_1|-d-\alpha \lambda }{\alpha +\Delta t^2(w_1+w_2)}$$ $$\Delta {\mathbf{p}}_1=+w_1\Delta t^2\Delta \lambda \mathbf{n}$$ $$\Delta {\mathbf{p}}_2=-w_2\Delta t^2\Delta \lambda \mathbf{n}$$

7.5 数值例子

参数：$m_1=m_2=1$，$\Delta t=0.02$，$\alpha =0.1$，$d=2$，${\mathbf{p}}_1=(0,0)$，${\mathbf{p}}_2=(3,0)$。

初始：$\lambda =0$，$C=1$，$\mathbf{n}=(1,0)$，$\text{denom}=2$。

$$\Delta \lambda =\frac{1-0}{0.1+0.0004\times 2}=\frac{1}{0.1008}\approx 9.92$$ $$\Delta {\mathbf{p}}_1=1\times 0.0004\times 9.92\times (1,0)=(0.00397,0)$$ $$\Delta {\mathbf{p}}_2=-1\times 0.0004\times 9.92\times (1,0)=(-0.00397,0)$$

新位置：${\mathbf{p}}_1=(0.00397,0)$，${\mathbf{p}}_2=(2.99603,0)$，新距离 $\approx 2.992$。

---

第八部分：完整算法总结

初始化: p, v, λ = 0

每一时间步 Δt:
    # 1. 预测位置
    p* = p + Δt·v + Δt²·M⁻¹·f_ext
    
    # 2. 迭代求解约束
    for iteration = 1 to solver_iterations:
        for each constraint j:
            # 计算当前值
            C = C_j(p)
            ∇C = gradient(C_j, p)
            denom = Σ w_i ||∇C_i||²
            
            # 更新 λ
            Δλ = (C - α_j·λ_j) / (α_j + Δt²·denom)
            λ_j = λ_j + Δλ
            
            # 更新位置
            for each particle i in constraint j:
                p_i = p_i - w_i·Δt²·Δλ·∇C_i
    
    # 3. 更新速度
    v = (p - p_old) / Δt

---

附录：概念速查表

符号	含义	公式
$\mathbf{p}$	位置	-
$\mathbf{v}$	速度	-
$\mathbf{M}$	质量矩阵	$\text{diag}(m_i)$
$w_i$	逆质量	$1/m_i$
$\Delta t$	时间步长	-
${\mathbf{p}}^{\ast }$	预测位置	${\mathbf{p}}^n+\Delta t{\mathbf{v}}^n$
$C(\mathbf{p})$	约束函数	如 $
$\alpha$	柔度	$1/k$（刚度倒数）
$\lambda$	拉格朗日乘子	约束力大小
$\text{denom}$	分母	$\sum w_i\Vert \nabla C{\Vert }^2$

---

文档结束