1. 从一个弹簧说起为什么积分方式决定了仿真的生死如果你做过布料模拟、头发摆动、软体形变这类效果大概率绕不开弹簧质点系统Mass-Spring System。它的思路极其朴素把物体拆成一堆带质量的点点与点之间用弹簧连接每个点受重力、弹簧力、阻尼力然后靠牛顿第二定律算出加速度再一步步推进时间得到下一帧的位置。听起来简单到不像话但真正上手写过的人都知道坑几乎全在“怎么推进时间”这一步上。同样一套弹簧质点模型用不同的积分方式去解结果可能是稳定的布料也可能是瞬间炸成一团乱麻的粒子云。这中间的差别就是显式积分和隐式积分的分水岭。这篇内容我打算把这两类积分方式从头推一遍公式不跳步然后把它们落到弹簧质点系统上最后给出一份可以直接跑的参考实现思路。适合正在做物理仿真、图形学入门、或者单纯想搞明白“为什么我的布料一跑就爆”的人。哪怕你数学忘得差不多了我也会尽量用生活化的方式把每一步讲清楚。先说结论方便你带着预期往下看显式积分实现简单、每步计算便宜但对时间步长极其敏感弹簧越硬越容易炸隐式积分每步计算贵、要解线性系统但可以用大得多的时间步长稳定性强得多。弹簧质点系统里弹簧刚度stiffness一高显式基本就没法用了这时候隐式几乎是唯一选择。下面我会先讲清楚弹簧质点系统的力学建模再分别推导显式和隐式积分接着讲隐式里最关键的线性系统怎么解最后给一份结构清晰的参考代码和一堆踩坑经验。2. 弹簧质点系统的力学建模先把“力”算明白2.1 质点、弹簧与受力分析弹簧质点系统的基本单元就两样东西质点有质量、有位置、有速度和弹簧连接两个质点有原长、有刚度、有阻尼。对任意一个质点 i它受到的力来自三部分重力F_g m_i * gg 通常取 (0, -9.8, 0)。弹簧力来自所有与它相连的弹簧。阻尼力抑制振荡让系统能量慢慢耗散。弹簧力怎么算设弹簧连接质点 i 和 j当前长度向量是d x_j - x_i l |d|弹簧原长记作 l_0刚度记作 k_s。那么弹簧的弹性力大小满足胡克定律f_s k_s * (l - l_0)方向沿着 d 的单位向量。注意符号如果 l l_0被拉长弹簧要往回拉所以作用在 i 上的力方向是 d 方向把 i 拉向 j作用在 j 上的力方向相反。写成向量F_i k_s * (l - l_0) * (d / l) F_j - k_s * (l - l_0) * (d / l)这里有个细节很多人第一次写会搞反当 l l_0 时(l - l_0) 为正i 受到的力应该指向 j也就是 d 方向所以 F_i 加的是正号。你可以拿一根被拉长的橡皮筋在脑子里过一遍确认方向没错。阻尼力通常用相对速度来算。设阻尼系数 k_d相对速度v_rel (v_j - v_i)阻尼力沿 d 方向的分量f_d k_d * dot(v_rel, d / l) F_i k_d * dot(v_rel, d / l) * (d / l) F_j - k_d * dot(v_rel, d / l) * (d / l)把所有弹簧的贡献累加再加上重力和全局阻尼就得到每个质点当前受到的合力 F_i。2.2 从力到加速度牛顿第二定律有了合力加速度就是a_i F_i / m_i然后我们要做的是已知当前时刻 t 的位置 x(t) 和速度 v(t)求下一时刻 tΔt 的位置和速度。这就是时间积分要解决的问题。连续时间下运动方程是dx/dt v dv/dt a F(x, v) / m注意 F 同时依赖位置弹簧力和速度阻尼力所以这是一个耦合的微分方程组。数值积分就是用离散的方式近似求解它。2.3 为什么“刚度”是稳定性的命门在往下讲积分之前必须先建立一个直觉弹簧越硬k_s 越大系统里最快的振荡频率就越高对时间步长的要求就越苛刻。单个弹簧质点振子的角频率大约是ω sqrt(k_s / m)周期 T 2π / ω。显式积分要稳定时间步长 Δt 必须小于某个和 T 相关的阈值后面会具体推。k_s 增大 10 倍ω 增大 sqrt(10) ≈ 3.16 倍允许的 Δt 就缩小到约 1/3.16。布料模拟里为了不让布被拉长k_s 往往设得很大于是显式积分的 Δt 被压到极小一帧要算几千上万步根本跑不动。这就是为什么工业级布料、软体仿真几乎都用隐式积分。理解了这一点后面的推导你才知道自己在解决什么问题。3. 显式积分简单直接但有个致命脾气3.1 显式欧拉最朴素的推进方式显式欧拉Explicit Euler是最容易理解的积分方式。它的逻辑是用当前时刻的力算出加速度然后假设这个加速度在 Δt 内不变直接往前推。公式v(tΔt) v(t) a(t) * Δt x(tΔt) x(t) v(t) * Δt注意这里位置更新用的是旧速度v(t)这是显式欧拉的一个标志。也有人用更新后的速度去更新位置v(tΔt) v(t) a(t) * Δt x(tΔt) x(t) v(tΔt) * Δt这叫半隐式欧拉Semi-implicit Euler也叫辛欧拉。别被名字骗了它本质上还是显式的因为 a(t) 只依赖当前状态。半隐式欧拉在物理仿真里比显式欧拉更常用因为它对振荡系统的能量表现更好稳定性也略强。3.2 显式积分的稳定性推导拿最简单的单弹簧振子来分析。设质点质量 m弹簧刚度 k忽略阻尼和重力位置偏离平衡点 x则a -(k/m) * x记 ω² k/m则 a -ω² x。用半隐式欧拉v_{n1} v_n - ω² x_n Δt x_{n1} x_n v_{n1} Δt把第一式代入第二式x_{n1} x_n (v_n - ω² x_n Δt) Δt x_n v_n Δt - ω² x_n Δt²再写出 v_{n1} 的表达式可以把它写成一个矩阵形式[ x_{n1} ] [ 1 - ω²Δt² Δt ] [ x_n ] [ v_{n1} ] [ -ω²Δt 1 ] [ v_n ]系统稳定要求这个矩阵的谱半径 ≤ 1。经过推导这里省略特征值求解的中间步骤结论是标准的稳定条件是Δt 2 / ω 2 * sqrt(m / k)这就是显式积分的CFL 类条件。翻译成人话弹簧越硬、质量越小允许的时间步长越小。举个具体数字。假设 m 1k 10000则 ω 100Δt 0.02 秒。如果 k 1000000ω 1000Δt 0.002 秒。布料模拟里 k 动辄几十万Δt 就得压到毫秒级以下一帧 16.7 毫秒要跑几十上百个子步。这就是显式积分在硬弹簧面前的真实处境。3.3 显式积分的实操要点与坑写显式积分时有几个地方特别容易翻车先更新速度再更新位置。用半隐式欧拉别用纯显式欧拉否则能量会莫名其妙增长系统更容易炸。子步substep是必须的。不要指望一个 Δt 走完一帧通常要把一帧切成 N 个子步每个子步重新算力。N 取多少按上面的稳定条件反推Δt_sub ≤ 2*sqrt(m/k)然后 N ceil(frame_dt / Δt_sub)。阻尼别设太大也别设太小。阻尼太小振荡不衰减看起来像果冻乱抖阻尼太大整个系统像泡在糖浆里动作发黏。一般从 k_d 0.1 * k_s 这个量级开始调。力的方向符号一定要验证。我见过太多人因为弹簧力方向写反导致质点互相排斥整个网格瞬间飞散。提示显式积分调试时先把重力关掉、只留一根弹簧观察两个质点能不能稳定振荡。这一步过了再逐步加复杂度。4. 隐式积分用“未来”的力来换稳定性4.1 隐式欧拉的核心思想显式积分用当前时刻的力隐式积分则用下一时刻的力v(tΔt) v(t) a(tΔt) * Δt x(tΔt) x(t) v(tΔt) * Δt问题来了a(tΔt) 依赖 x(tΔt)而 x(tΔt) 又依赖 v(tΔt)v(tΔt) 又依赖 a(tΔt)。这是一个隐式方程不能直接算出来必须解方程。为什么这样更稳定直觉上隐式积分相当于“预判了未来的力会把你拉回来”所以它对大步长更宽容。数学上隐式欧拉对线性系统是无条件稳定的A-stable不管 Δt 多大都不会发散。代价是每步要解一个方程组。4.2 一阶泰勒展开与线性化对 a(tΔt) 在 t 处做一阶泰勒展开a(tΔt) ≈ a(t) (∂a/∂x) * Δx (∂a/∂v) * Δv其中 Δx x(tΔt) - x(t)Δv v(tΔt) - v(t)。而由隐式欧拉Δx v(tΔt) * Δt (v(t) Δv) * Δt Δv a(tΔt) * Δt把 Δx 代入 a 的展开式a(tΔt) ≈ a(t) (∂a/∂x) * (v(t) Δv) * Δt (∂a/∂v) * Δv再代入 Δv a(tΔt) * Δt得到关于 a(tΔt) 的方程a(tΔt) a(t) (∂a/∂x) * v(t) * Δt (∂a/∂x) * a(tΔt) * Δt² (∂a/∂v) * a(tΔt) * Δt整理把含 a(tΔt) 的项移到一边a(tΔt) * (I - (∂a/∂x) * Δt² - (∂a/∂v) * Δt) a(t) (∂a/∂x) * v(t) * Δt这就是隐式积分的核心方程。左边括号里是一个矩阵右边是已知量。解出 a(tΔt)就能更新速度和位置。4.3 弹簧质点系统里的雅可比矩阵在弹簧质点系统里∂a/∂x 和 ∂a/∂v 是有明确物理意义的∂a/∂x力的位置导数除以质量反映弹簧刚度对位置的敏感度。对一根连接 i、j 的弹簧它贡献的块是一个 3x3 的子矩阵在三维里。∂a/∂v力的速度导数除以质量反映阻尼对速度的敏感度。对弹簧力 F k_s * (l - l_0) * (d/l)它对位置的导数推导起来有点繁琐但结论是∂F/∂x_j k_s * (I - (l_0/l) * (I - d⊗d/l²))其中 ⊗ 是外积I 是单位矩阵。这个矩阵在 l ≈ l_0 时接近 k_s * I在 l 远大于 l_0 时会变小。阻尼力的导数更简单大致是 k_d * (d⊗d/l²)。把所有弹簧的贡献组装起来就得到一个大的稀疏矩阵。这个矩阵通常记作A I - (∂a/∂x) * Δt² - (∂a/∂v) * Δt然后解线性系统A * a(tΔt) b其中 b a(t) (∂a/∂x) * v(t) * Δt。4.4 为什么隐式积分“贵但值”隐式积分每步要解一个 N×N 的线性系统N 是质点数的 3 倍因为每个质点三个分量。直接求逆是 O(N³)对几千个质点根本不可行。所以实际都用迭代法比如共轭梯度CG或者预条件共轭梯度PCG利用矩阵的稀疏性每步复杂度接近 O(N)。即便如此隐式每步的计算量还是比显式大不少。但关键在于隐式可以用大得多的 Δt。显式可能要 100 个子步隐式 1 到 5 个子步就够了。算总账隐式反而更快而且稳定。注意隐式积分不是“绝对不会炸”。如果线性系统解得不收敛或者雅可比矩阵近似得太粗糙照样会出问题。稳定性是有条件的只是条件比显式宽松得多。5. 线性系统求解隐式积分的真正战场5.1 共轭梯度法为什么适合这里隐式积分要解的矩阵 A 有几个好性质稀疏每个质点只和相邻质点有连接矩阵非零元很少。对称正定在合理建模下弹簧力是保守力阻尼也通常建模成对称的。对称正定加稀疏简直是共轭梯度法CG的完美猎物。CG 不需要存储整个矩阵只需要能计算“矩阵乘向量”这个操作。对弹簧质点系统来说矩阵乘向量可以逐弹簧、逐质点地算内存开销极小。CG 的迭代公式这里不展开全部推导核心是它在一个 Krylov 子空间里找最优解理论上 N 步收敛N 是矩阵维度但实际因为条件数好往往几十步就够。5.2 预条件让 CG 快起来的秘诀CG 的收敛速度取决于矩阵的条件数。弹簧质点系统的矩阵条件数往往不小硬弹簧和软弹簧混在一起所以需要预条件。最常用的预条件子是对角预条件Jacobi 预条件取 A 的对角线元素构造一个对角矩阵 M用 M⁻¹ 去“缩放”问题。实现简单效果通常也不错。更高级的还有块对角预条件、不完全 Cholesky 分解等但实现复杂度上升。我的经验是先上 Jacobi 预条件如果收敛太慢再考虑更复杂的方案。5.3 迭代终止条件怎么定CG 迭代不能无限跑得有个终止条件。常用的是相对残差||r|| / ||b|| toltol 一般取 1e-4 到 1e-6。tol 太小浪费时间tol 太大解不准系统可能抖动。还有一个坑每帧的线性系统是变化的但变化不大。可以用上一帧的解作为这一帧的初始猜测warm start能显著减少迭代次数。这个技巧在实时仿真里几乎是标配。6. 参考实现一份能跑通的弹簧质点仿真骨架6.1 数据结构设计先定义质点class Particle: def __init__(self, pos, mass, pinnedFalse): self.pos np.array(pos, dtypefloat) self.vel np.zeros(3) self.force np.zeros(3) self.mass mass self.inv_mass 0.0 if pinned else 1.0 / mass self.pinned pinned弹簧class Spring: def __init__(self, i, j, rest_len, ks, kd): self.i i self.j j self.rest_len rest_len self.ks ks self.kd kd用 inv_mass 而不是 mass是为了方便处理固定点pinned固定点 inv_mass 0力再大也不动。6.2 显式积分主循环def step_explicit(particles, springs, dt, gravity): # 1. 清零力加重力 for p in particles: p.force gravity * p.mass if not p.pinned else np.zeros(3) # 2. 累加弹簧力 for s in springs: pi, pj particles[s.i], particles[s.j] d pj.pos - pi.pos l np.linalg.norm(d) if l 1e-12: continue dir d / l # 弹性力 fs s.ks * (l - s.rest_len) # 阻尼力 v_rel pj.vel - pi.vel fd s.kd * np.dot(v_rel, dir) f (fs fd) * dir pi.force f pj.force - f # 3. 半隐式欧拉更新 for p in particles: if p.pinned: continue acc p.force * p.inv_mass p.vel acc * dt p.pos p.vel * dt这段代码逻辑清晰但注意它只在 dt 足够小时稳定。实际用的时候外面要套一层子步循环。6.3 隐式积分主循环简化版隐式版本复杂得多这里给一个结构骨架重点展示线性系统怎么组装和求解def step_implicit(particles, springs, dt, gravity, cg_iters50): n len(particles) # 1. 计算当前加速度 a(t) compute_forces(particles, springs, gravity) a0 np.array([p.force * p.inv_mass for p in particles]) # 2. 组装右端项 b a0 (da/dx) * v * dt # 这里 (da/dx)*v 需要逐弹簧计算略去细节 b a0 compute_dadx_times_v(particles, springs) * dt # 3. 组装矩阵 A 的对角线用于 Jacobi 预条件 diag compute_matrix_diag(particles, springs, dt) # 4. 用 CG 解 A * a_new b a_new conjugate_gradient( lambda x: matvec_A(x, particles, springs, dt), b, diag, iterscg_iters ) # 5. 更新速度和位置 for idx, p in enumerate(particles): if p.pinned: continue p.vel a_new[idx] * dt p.pos p.vel * dt其中 matvec_A 是矩阵乘向量操作逐弹簧累加贡献。这个函数是性能热点要尽量写高效。6.4 参数选择与调优经验参数显式推荐隐式推荐说明时间步长 Δt1e-4 ~ 1e-3 秒1e-2 ~ 1.6e-2 秒隐式可用整帧步长子步数20 ~ 2001 ~ 5显式必须多子步弹簧刚度 k_s1e3 ~ 1e41e5 ~ 1e7隐式能扛硬弹簧阻尼 k_d0.1 * k_s0.01 * k_s隐式阻尼可小些CG 迭代数不适用20 ~ 100配合 warm startCG 容差不适用1e-4太严浪费时间这张表是我自己调参时总结的具体数值要看场景。布料模拟 k_s 通常要 1e5 以上显式基本没戏。7. 常见问题与排查技巧实录7.1 系统爆炸的几种典型原因现象一一帧之内质点飞到无穷远。最常见的原因是时间步长超过稳定极限。排查方法把 k_s 调小 10 倍如果稳定了就是步长问题。解决减小子步 Δt或者换隐式。现象二系统缓慢膨胀像充气。这通常是显式欧拉非半隐式的能量增长问题。解决改用半隐式欧拉先更新速度再更新位置。现象三网格抖动、高频振荡。阻尼不足。解决增大 k_d或者加全局速度阻尼。现象四隐式积分也炸。检查 CG 是否收敛。如果 CG 没收敛就用了错误的解系统会发散。解决增大迭代数或者放宽容差或者检查矩阵组装是否有 bug。7.2 性能优化清单用空间局部性质点位置、速度用连续数组存储SoA 布局比对象数组快得多。避免重复开方l sqrt(dot(d,d))能复用就复用。CG 的 warm start用上一帧的解做初值迭代数能减半。固定点提前剔除pinned 质点不参与线性系统能减小矩阵规模。并行化弹簧力累加和矩阵乘向量都可以并行但要注意写冲突用原子操作或分块。7.3 一个容易被忽略的坑弹簧原长的初始化很多人建完网格直接设 rest_len 当前长度这没错。但如果初始网格本身就有拉伸比如从别的格式导入的模型rest_len 设成当前长度会导致初始状态就有应力系统一开始就抖。正确做法是要么确保初始网格无拉伸要么把 rest_len 设成设计长度而不是当前长度。8. 显式与隐式的选择一张决策表场景推荐方案理由弹簧刚度低、实时性要求高显式 多子步实现简单够用弹簧刚度高布料、软体隐式 CG显式步长太小跑不动质点数量少 100显式隐式的线性系统开销不划算质点数量多 1000隐式显式子步数会爆炸需要精确能量守恒隐式或辛积分器显式欧拉能量不守恒快速原型验证显式代码量小调试快我个人的经验是先用显式把逻辑跑通确认物理建模没问题再换隐式解决稳定性。一上来就写隐式调试难度会大很多因为你分不清是建模错了还是线性系统解错了。9. 从公式到代码几个推导细节的补充9.1 弹簧力雅可比矩阵的完整推导前面给了 ∂F/∂x_j 的结论这里补一下推导思路。弹簧力F k_s * (l - l_0) * (d / l)其中 d x_j - x_il |d|。对 x_j 求导要用到链式法则。关键是∂l/∂d d / l ∂(d/l)/∂d (I - d⊗d/l²) / l代入后整理得到∂F/∂x_j k_s * (d⊗d/l²) k_s * (l - l_0)/l * (I - d⊗d/l²)化简∂F/∂x_j k_s * (I - (l_0/l) * (I - d⊗d/l²))这个矩阵是对称的符合物理直觉作用力反作用力。在 l l_0 时它等于 k_s * I就是纯刚度。9.2 为什么用加速度而不是位置作为未知量隐式积分里未知量可以选位置、速度或加速度。选加速度的好处是方程形式最简洁右端项直接是当前加速度矩阵 A 的形式也规整。选位置的话方程里会出现 Δt² 的系数数值上更容易出问题。9.3 阻尼的隐式处理阻尼力依赖速度在隐式框架下速度也是未知的。所以阻尼的雅可比 ∂a/∂v 会进入矩阵 A。这也是为什么 A 里有一项 -(∂a/∂v) * Δt。如果忽略这一项把阻尼显式处理稳定性会打折扣但实现简单。实际中很多人这么做效果也能接受。10. 我踩过的坑和最后的小技巧写弹簧质点仿真这些年有几个教训是文档里不会写的。第一别迷信“无条件稳定”。隐式积分的无条件稳定是对线性系统说的弹簧质点系统是非线性的大步长下线性化误差会累积照样可能出问题。我的做法是隐式也用适中的步长别一上来就整帧步长。第二CG 的收敛性比想象中脆弱。如果弹簧刚度差异极大比如一根 k1e7一根 k1e2矩阵条件数会很差CG 收敛慢甚至不收敛。这时候要么用更强的预条件要么把刚度差异控制住。第三调试时把可视化做扎实。光看数字很难判断系统状态。把质点画出来用颜色表示速度大小一眼就能看出哪里在炸。第四warm start 是免费的加速。上一帧的解作为这一帧 CG 的初值几乎不增加代码复杂度但迭代数能降不少。这个技巧我强烈建议默认开启。最后分享一个小技巧如果你在显式和隐式之间纠结可以先写一个统一的框架把积分方式做成可切换的模块。这样调试时能快速对比两种方式的表现也方便后续扩展。我自己的代码里就是这么做的切换积分器只需要改一个枚举值。这个方向后续还能扩展的地方很多比如加入弯曲弹簧bending spring让布料更挺括或者用位置约束position-based dynamics替代力-based 方法稳定性又是另一套逻辑。但不管怎么扩展把显式和隐式这两套积分的地基打牢后面都是水到渠成的事。
