1. 简介
布料仿真中,我们通常将布料剖分为三角形网格(或四边形网格),并用弹簧-质点模型构造动力学系统:质点即三角形的顶点,弹簧即三角形的边。质点在外力(如,重力)和内力(弹簧力)的作用下根据牛顿第二定律运动:
{x˙=vMv˙=Fexternal+Finternal
这就是我们需要求解的常微分方程,通常我们很难求其解析解,只能求其数值解。常用方法有:显式/隐式欧拉法,Symplectic Euler,Midpoint method,Leapfrog integration等。
2. 显式欧拉法(Explicit Euler)
形如,
{xn+1=xn+Δtvnvn+1=vn+M−1Fn
称为显式积分,可以看出,位置和速度的积分用的都是显式方法。显式积分需控制积分步长 Δt ,以保证数值求解的稳定性。 这里例证了显式条件稳定性。
3. 隐式欧拉法(Implicit Euler)
形如,
{xn+1=xn+Δtvn+1vn+1=vn+M−1Fn+1
称为显式积分,可以看出,位置和速度的积分用的都是隐式方法。 从理论上说,隐式欧拉法是无条件稳定的。即,积分步长可以任意大,数值求解都稳定(稳定不意味着精确)。然而, 隐式方法需要求解非线性方程,通常要线性近似,因此实际中,隐式方法并不是无条件稳定的。即,仍要控制积分步长以保证数值求解的稳定。
这里例证了隐式欧拉法的无条件稳定性。
另注,隐式欧拉法并不是唯一的隐式方法,例如下面将要介绍的midpoint method也属于隐式方法。只是由于其简单及良好性能,隐式欧拉应用得更广泛。
4. Symplectic Euler
形如(1),
{xn+1=xn+Δtvnvn+1=vn+M−1Fn+1
或(2),
{xn+1=xn+Δtvn+1vn+1=vn+M−1Fn
称为Symplectic Euler。
- 可以看出(1)位置的更新使用显式方法,速度的更新使用隐式方法;(2)位置的更新使用隐式方法,速度的更新使用显式方法。因此Symplectic Euler又称为semi-implicit Euler。
- 由于(1)与隐式欧拉法的计算量相同,因此实际中很少使用(1),而直接使用隐式欧拉法;
- (2)与显式欧拉法的计算量相同,但比显式欧拉法更稳定,可以使用更大的积分步长,因此实际中很少使用显式欧拉法,而使用(2)。从理论上分析,仿真谐振动系统(harmonic oscillation system),只要保证 Δt<2ω=2m/k−−−−√ ,形式(2)的Symplectic Euler就是稳定的。
5. 中点法(Midpoint method)
形如,
yn+1=yn+hf(tn+h2,yn+h2f(tn,yn))
称为 显式中点法(Explicit Midpoint method),又称为改进的欧拉法(modified Euler method)。
形如,
yn+1=yn+hf(tn+h2,12(yn+yn+1))
称为 隐式中点法(Implicit Midpoint method)。
由上面两个式子可以看出, f 总是计算在
由两种方法 估计 f(tn+h2,y(tn+h2)) 的值:
1. 将 y(tn+h2) 在 tn 处做一阶泰勒展开,得,
y(tn+h2)≈y(tn)+h2f(tn,y(tn))=yn+h2f(tn,yn)
2. 用 y(tn) 和 y(tn+h) 的均值近似:
y(tn+h2)≈12(y(tn)+y(tn+h))=12(yn+yn+1)
- 中点法的局部误差为 O(h3) ,全局误差为 O(h2) 。因此比欧拉法的 O(h2) 局部误差, O(h) 全局误差要精确。当然,隐式欧拉法的计算量更大。
- 中点法实际上是Runge-Kutta method(以后会专门介绍)的一种特例。
以上几种方法又可以写成如下通用形式:
{1hΔx=(1−β)vn+βvn+11hMΔv=(1−α)Fn+αFn+1
- 当 α=β=0 时,为显式欧拉法
- 当 α=β=1 时,为隐式欧拉法
- 当 α=0,β=1 或 α=1,β=0 时,为Symplectic Euler(Semi-implicit Euler)。
- 当 α=β=12 时,参考论文3 Cloth simulation using soft constraints 错误的将其认为是隐式中点法,而实际上中点法的核心是求 tn+h2 时刻的导数 Fn+12 ,并不是 tn 时刻和 tn+1 时刻导数的均值 12(Fn+Fn+1) 。当 α=β=12 时,应称为 显式隐式混合法(mixed explicit-implicit),实际上除前三种情况外,都可以称为显式隐式混合法。
6. Leapfrog method
形如,
⎧⎩⎨xn+1=xn+Δtvn+12vn+12=vn−12+M−1ΔtFn
称为Leapfrog method。此式又可改写成如下形式:
{xn+1=xn+Δtvn+12M−1FnΔt2vn+1=vn+12M−1(Fn+Fn+1)Δt
由此可知,Leapfrog method也是一种隐式方法。另外,Leapfrog method的局部误差为 O(h3) ,全局误差为 O(h2) 。只要时间步长 Δt≤1/ω ,该方法就是稳定的。
7. 总结
- 一般来说显式法稳定性差,需要较小的积分步长,但能更好的保持能量(energy preserving);
- 隐式法稳定性好,可以使用更大的积分步长,但会增加能量耗散(energy dissipation)。
- 欧拉法(显式和隐式),Symplectic Euler的局部误差为 O(h2) ,全局误差为 O(h) ;中点法,Leapfrog method的局部误差为 O(h3) ,全局误差为 O(h2) 。
- 有时候我们会对各种力分别做显示和隐式处理,称为IMEX方法。
- 一般情况下,减小时间步长 Δt ,显式欧拉误差的下降速度要稍快与隐式欧拉。
8. 参考资料
- 维基百科 Leapfrog integration词条,Semi implicit Euler method词条,Midpoint method词条。
- Interactive simulation of elastic deformable materials
Servin M, Lacoursiere C, Melin N. SIGRAD 2006.- Frâncu M, Moldoveanu F. Cloth simulation using soft constraints[J]. 2015.