企业官网建设流程全解析

1. 从一个弹簧说起：为什么积分方式决定了仿真的生死

如果你做过布料模拟、头发摆动、软体形变这类效果，大概率绕不开弹簧质点系统（Mass-Spring System）。它的思路极其朴素：把物体拆成一堆带质量的点，点与点之间用弹簧连接，每个点受重力、弹簧力、阻尼力，然后靠牛顿第二定律算出加速度，再一步步推进时间，得到下一帧的位置。

听起来简单到不像话，但真正上手写过的人都知道，坑几乎全在“怎么推进时间”这一步上。同样一套弹簧质点模型，用不同的积分方式去解，结果可能是稳定的布料，也可能是瞬间炸成一团乱麻的粒子云。这中间的差别，就是显式积分和隐式积分的分水岭。

这篇内容我打算把这两类积分方式从头推一遍，公式不跳步，然后把它们落到弹簧质点系统上，最后给出一份可以直接跑的参考实现思路。适合正在做物理仿真、图形学入门、或者单纯想搞明白“为什么我的布料一跑就爆”的人。哪怕你数学忘得差不多了，我也会尽量用生活化的方式把每一步讲清楚。

先说结论，方便你带着预期往下看：显式积分实现简单、每步计算便宜，但对时间步长极其敏感，弹簧越硬越容易炸；隐式积分每步计算贵、要解线性系统，但可以用大得多的时间步长，稳定性强得多。弹簧质点系统里，弹簧刚度（stiffness）一高，显式基本就没法用了，这时候隐式几乎是唯一选择。

下面我会先讲清楚弹簧质点系统的力学建模，再分别推导显式和隐式积分，接着讲隐式里最关键的线性系统怎么解，最后给一份结构清晰的参考代码和一堆踩坑经验。

2. 弹簧质点系统的力学建模：先把“力”算明白

2.1 质点、弹簧与受力分析

弹簧质点系统的基本单元就两样东西：质点（有质量、有位置、有速度）和弹簧（连接两个质点，有原长、有刚度、有阻尼）。

对任意一个质点 i，它受到的力来自三部分：

• 重力：F_g = m_i * g，g 通常取 (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_{n+1} = v_n - ω² x_n Δt x_{n+1} = x_n + v_{n+1} Δt

把第一式代入第二式：

x_{n+1} = x_n + (v_n - ω² x_n Δt) Δt = x_n + v_n Δt - ω² x_n Δt²

再写出 v_{n+1} 的表达式，可以把它写成一个矩阵形式：

[ x_{n+1} ] [ 1 - ω²Δt² Δt ] [ x_n ] [ v_{n+1} ] = [ -ω²Δt 1 ] [ v_n ]

系统稳定要求这个矩阵的谱半径 ≤ 1。经过推导（这里省略特征值求解的中间步骤，结论是标准的），稳定条件是：

Δt < 2 / ω = 2 * sqrt(m / k)

这就是显式积分的CFL 类条件。翻译成人话：弹簧越硬、质量越小，允许的时间步长越小。

举个具体数字。假设 m = 1，k = 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|| < tol

tol 一般取 1e-4 到 1e-6。tol 太小，浪费时间；tol 太大，解不准，系统可能抖动。

还有一个坑：每帧的线性系统是变化的，但变化不大。可以用上一帧的解作为这一帧的初始猜测（warm start），能显著减少迭代次数。这个技巧在实时仿真里几乎是标配。

6. 参考实现：一份能跑通的弹簧质点仿真骨架

6.1 数据结构设计

先定义质点：

class Particle: def __init__(self, pos, mass, pinned=False): self.pos = np.array(pos, dtype=float) 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_iters=50): 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, iters=cg_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 参数选择与调优经验

这张表是我自己调参时总结的，具体数值要看场景。布料模拟 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. 显式与隐式的选择：一张决策表

我个人的经验是：先用显式把逻辑跑通，确认物理建模没问题，再换隐式解决稳定性。一上来就写隐式，调试难度会大很多，因为你分不清是建模错了还是线性系统解错了。

9. 从公式到代码：几个推导细节的补充

9.1 弹簧力雅可比矩阵的完整推导

前面给了 ∂F/∂x_j 的结论，这里补一下推导思路。弹簧力：

F = k_s * (l - l_0) * (d / l)

其中 d = x_j - x_i，l = |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 的收敛性比想象中脆弱。如果弹簧刚度差异极大（比如一根 k=1e7，一根 k=1e2），矩阵条件数会很差，CG 收敛慢甚至不收敛。这时候要么用更强的预条件，要么把刚度差异控制住。

第三，调试时把可视化做扎实。光看数字很难判断系统状态。把质点画出来，用颜色表示速度大小，一眼就能看出哪里在炸。

第四，warm start 是免费的加速。上一帧的解作为这一帧 CG 的初值，几乎不增加代码复杂度，但迭代数能降不少。这个技巧我强烈建议默认开启。

最后分享一个小技巧：如果你在显式和隐式之间纠结，可以先写一个统一的框架，把积分方式做成可切换的模块。这样调试时能快速对比两种方式的表现，也方便后续扩展。我自己的代码里就是这么做的，切换积分器只需要改一个枚举值。

这个方向后续还能扩展的地方很多，比如加入弯曲弹簧（bending spring）让布料更挺括，或者用位置约束（position-based dynamics）替代力-based 方法，稳定性又是另一套逻辑。但不管怎么扩展，把显式和隐式这两套积分的地基打牢，后面都是水到渠成的事。