一个无阻尼弹簧从拉伸位置释放，理论上应一直以相同振幅振荡。用最直接的 Euler 更新位置和速度，程序却越振越大；把位置改成使用更新后的速度，画面稳定了，但它是否真的保存了能量？

本篇固定一个能写出解析解的单弹簧，比较两种更新顺序。所有实验到达同一个 11 秒时刻，同时观察位移、速度、能量和步长。半隐式 Euler 只在特定步长范围内有界，能量有界也不代表运动相位准确。

动画姿态与动力学状态有不同输入

第 29 篇用关键帧和绝对时间求姿态，时间一给定就能插值得到结果。动力学模型则保存状态，下一状态依赖当前状态、受力和积分方法。只知道位置，通常还不够决定下一步向哪里运动。

取一维质点，弹簧固定端与静止平衡点位于 x=0。x 是偏离平衡点的有符号位移，v 是速度；当前例子没有重力、阻尼、碰撞和弹簧长度下限。它不是一条只能承受拉力的绳子，x<0 时仍按同一 Hooke 定律回复。

量纲具体规定为米、千克、秒。质量 m=1 kg，弹簧刚度 k=4 N/m，受力 F=−kx，正位移产生负向力，负位移产生正向力。由 F=ma 得到：

x˙=v,v˙=−kmx=−ω2x,ω=k/m=2 rad/s.\dot x=v,\qquad \dot v=-\frac{k}{m}x=-\omega^2x,\qquad \omega=\sqrt{k/m}=2\ {\rm rad/s}.

点号表示对时间求导。速度是位置变化率，加速度是速度变化率；x 与 v 的单位不同，不能把它们当作两项同单位坐标直接比较欧氏误差。下面分别报告位移误差（米）和速度误差（米/秒）。

初态 x(0)=1、v(0)=0，没有随机扰动。独立解析解与真实机械能为：

x(t)=cos⁡(2t),v(t)=−2sin⁡(2t),E=12mv2+12kx2=2 J.x(t)=\cos(2t),\qquad v(t)=-2\sin(2t),\qquad E=\tfrac12mv^2+\tfrac12kx^2=2\ {\rm J}.

能量恒定可以由求导检查：dE/dt=mv·v̇+kx·ẋ=−kxv+kxv=0。这是连续微分方程的性质，还没有说明离散算法也会保存它。

Euler 的一个步长近似了什么

若时间步长为 h，Taylor 展开给出 x(t+h)=x(t)+h·ẋ(t)+O(h²)。显式 Euler 保留前两项，位置和速度都使用步长起点的信息：

xn+1=xn+hvn,vn+1=vn−hω2xn.x_{n+1}=x_n+hv_n,\qquad v_{n+1}=v_n-h\omega^2x_n.

局部截断误差是 O(h²)，在固定有限时长与适当光滑条件下，全局误差是一阶。这个精度说法不保证长期稳定；在一个振荡系统里，每步很小的结构性误差仍可能持续积累。

当前半隐式 Euler 选择先更新速度，再用新速度推进位置：

vn+1=vn−hω2xn,xn+1=xn+hvn+1.v_{n+1}=v_n-h\omega^2x_n,\qquad x_{n+1}=x_n+hv_{n+1}.

这也是 symplectic Euler 的一种次序。另有先位置后速度的变体，公式和对应修改能量不同，不能只写“半隐式”而不说明执行顺序。当前力只依赖旧位置，速度更新仍然直接计算，没有求解一个隐式方程组，也没有使用完整隐式 Euler。

以 h=0.05 秒走第一步，显式得到 (x,v)=(1,−0.2)，半隐式得到 (0.99,−0.2)。如果 C++ 原地赋值后又无意读到了新变量，可能把两种方法混成另一套算法；spring31.hpp 按值接收状态，明确保留这两种顺序。

以下矩阵与离散解代入固定的 SI 参数：h 取步长的秒数，x 取米数，v 取米/秒数，能量表达式的结果为焦耳数。一般模型的步长条件仍用无量纲 hω 表示。

显式方法为什么一直增加能量

对本例，显式状态更新的矩阵为：

[xn+1vn+1]=[1h−4h1][xnvn].\begin{bmatrix}x_{n+1}\\v_{n+1}\end{bmatrix} =\begin{bmatrix}1&h\\-4h&1\end{bmatrix} \begin{bmatrix}x_n\\v_n\end{bmatrix}.

把两个新状态代入 E=2x²+v²/2，交叉项相互抵消，剩余项给出精确的离散增长律：

En+1=(1+4h2)En,En=2(1+4h2)n.E_{n+1}=(1+4h^2)E_n,\qquad E_n=2(1+4h^2)^n.

任意 h>0，因子都大于 1。因此当前无阻尼振子中，显式 Euler 的长期能量必然增长。减小 h 可以在有限时长内减轻误差，却没有将这个因子变成 1；不能将“小步长下暂时像振荡”解释为长期稳定。

程序对每一步的能量检查这个独立表达式。四个步长中最大相对差约 2.79×10⁻¹⁵。以 h=0.05 秒走 220 步，最终能量已从 2 J 增至 17.8538642 J，位移约 −2.98167 m，明显超出解析振幅 1 m。

这也解释了为什么不能为通过检查临时加阻尼：阻尼会更改微分方程，新的能量衰减不再是当前无阻尼算法的证据。当前实现保留不稳定结果，没有钳制位移或隐藏增大的能量。

半隐式的步长条件来自矩阵

将速度更新代入位置，半隐式矩阵为：

A=[1−4h2h−4h1],det⁡A=1,tr⁡A=2−4h2.A=\begin{bmatrix}1-4h^2&h\\-4h&1\end{bmatrix},\qquad \det A=1,\quad \operatorname{tr}A=2-4h^2.

当 0<h<1，两个特征值为一对不同的共轭复数，模长均为 1，因此该线性系统的离散轨迹有界。一般单弹簧的条件写成 0<hω<2；本例 ω=2，正好得到 h<1 秒。稳定性条件由当前矩阵推导，不从别的积分器套用一个同名阈值。

到 h=1 时，特征值重复为 −1，但矩阵不是 −I，不能只看模长就断言有界。对本篇初态，独立解为：

xn=(−1)n(1+2n),vn=(−1)n4n.x_n=(-1)^n(1+2n),\qquad v_n=(-1)^n4n.

振幅随 n 线性增长，能量随 n²增长。第 11 步得到 (−23,−44)，E=2026 J。这个临界点必须排除，不能把稳定条件写成 hω≤2。特殊初态可能没有同样增长，但方法的稳定性需要覆盖所有允许初态。

h>1 时，特征值成为互为倒数的实数，其中一个模长大于 1。h=1.1 秒、10 步后的能量约 4.14364×10⁸ J。无论采用哪一种更新顺序，都不应把“半隐式”理解为无需考虑刚度、质量与步长。

真实轨迹的四步长能量图，显式与半隐式在统一对数焦耳标度上比较

图中红线为显式，蓝线为半隐式，灰线为解析 E=2 J。纵轴是 log10(E/1 J)：0 对应 1 J，9 对应 10⁹ J；四面板使用同一标度。曲线连接实际离散采样，斜率不是连续瞬时功率，也不是线性能量变化率。

保存了一个二次量，不等于保存真实能量

对当前半隐式更新，直接展开还可以证明另一个量保持不变。带参数的形式是 Ẽ=E−hkxv/2，代入本例得到：

E~=12v2+2x2−2hxv.\widetilde E=\tfrac12v^2+2x^2-2hxv.

它包含位置与速度的交叉项，与真实 E 有区别。令无量纲 η=hω/2，本例 η 就是 h 的秒数。由 |hkxv/2|≤ηE，当 0<η<1 时有 (1−η)E≤Ẽ≤(1+η)E。初态的 Ẽ=2 J，从而得到能量界：

2 J1+η≤En≤2 J1−η.\frac{2\ {\rm J}}{1+\eta}\leq E_n\leq\frac{2\ {\rm J}}{1-\eta}.

这是当前模型与初态的界，不是所有物理系统的能量保证。h 达到 1 时下侧正定约束消失；h>1 时，即使 Ẽ仍在数学上守恒，也允许状态沿不定二次型方向越来越大。因此，“某个不变量没变”本身不能证明状态稳定。

h=0.05 的真实能量实际落在 1.90476193 到 2.10524812 J 之间，并非恒定 2 J。修改量最大绝对漂移约 2.89×10⁻¹⁵ J。h=1.1 的二次项已经很大，相减后的修改量绝对漂移约 5.00×10⁻⁸ J；按 max(2 J,E) 缩放后的误差约 1.21×10⁻¹⁶。记录两种误差，避免把大项消减产生的残差与稳定性混为一谈。

能量正确的时刻，位置仍可能错误

所有实验固定终点 T=11 秒，分别用 220、22、11、10 个完整步长。没有用不同总时间的末端结果比较精度，也没有为了对齐终点追加一个未说明的短步。

h，秒 半隐式终点 (x，v) 全程真实能量范围，J
0.05 (−0.998935，0.0361000) 1.90476–2.10525
0.5 (0，2) 2–4
1 (−23，−44) 2–2026
1.1 (12113.9，15548.0) 2–4.14364×10⁸

解析终点为 (−0.9999608264，0.0177026186)，E=2 J。h=0.5 的离散终点能量恰好也是 2 J，但位置却在 0 m，速度为 2 m/s：系统走到了另一相位。仅检验终点能量，甚至只检验有界，都无法保证运动准确。

稳定区的离散相位每步为 φ=2asin(h)，解析相位每步为 2h。对初态 (1,0)，离散解可写成：

xn=cos⁡(nϕ)−h1−h2sin⁡(nϕ),vn=−21−h2sin⁡(nϕ).x_n=\cos(n\phi)-\frac{h}{\sqrt{1-h^2}}\sin(n\phi),\qquad v_n=-\frac{2}{\sqrt{1-h^2}}\sin(n\phi).

当 h 很小时 φ接近 2h；h=0.5 时 φ=π/3，与解析每步 1 弧度已有差异。这个表达式还包含振幅与位置相位偏移，不能把所有误差都解释成频率变化。程序逐步与离散闭式比对，再与连续解析轨迹比较。

实际半隐式位移采样与解析振荡的同尺度对照，0.05秒和0.5秒步长

位移图使用同一米标度，灰色解析曲线取 441 个固定时刻，蓝色是实际半隐式采样折线。h=0.05 的全程最大位置差为 0.0585676 m，虽末端误差只有 0.00102583 m，全程不能按末端数值概括。连接粗采样的线段也不代表程序真的在两个状态之间积分出了线性运动。

验证与练习

完整 spring31.csv 与 spring31.json 保存八组轨迹的 534 条记录，包含每步状态、解析状态、真实/修改能量及各组误差。程序固定步数结束，没有随机数、性能排名或 GPU 代跑。

检查不仅覆盖初态静止，还用不同质量和非零初速度检查 k/m 与更新顺序。12 项非法输入拒绝非正质量/刚度/步长、非有限参数与状态；4 项溢出检查拒绝非有限状态或能量。首次检查曾主观要求小步长位置误差小于 0.04 m，实际值超出；改用离散闭式及推导误差界约 0.05924 m 后通过，没有更改积分算法或降低实验时间。

练习一：m=4 kg、k=4 N/m，当前速度先更新的半隐式 Euler 稳定步长范围是多少？h=2 秒能否包含在范围中？

解题要点：ω=1 rad/s，因此 0<h<2 秒。等号仍是重复特征根边界，不能包含。增大质量或减小刚度会扩大这一个线性振子的稳定步长上限，但不保证其他力或碰撞也稳定。

练习二：h=0.5 秒时，第 22 步 E=2 J，与解析完全一样。只保存一列能量日志是否足以检验动画轨迹？至少还应保存什么？

解题要点：不足。离散状态 (0,2) 与解析约 (−1,0.0177) 不同，应保存带时间戳的位移与速度并与相同时间的参考比较；全程误差和能量范围也应保留，避免终点巧合。

练习三：半隐式保持 Ẽ=2，为什么 h=1 时仍可以增长？代入 xₙ、vₙ验证第 11 步。

解题要点：h=1 时 Ẽ=0.5(v−2x)²，不再限制所有方向。x=−23、v=−44，v−2x=2，Ẽ=2；真实 E=0.5×44²+2×23²=2026 J。一个退化二次量无法约束真实能量。

1
2
make -C examples/computer-graphics build/spring31_check
examples/computer-graphics/build/spring31_check examples/computer-graphics/build

SpringState 只保存 x、v；积分函数按值返回下一状态，不隐藏累计时间。下一篇把固定模拟步长与显示帧率分开，并为碰撞说明步内发生的事件；改变显示频率不会自动改变这里的稳定条件。

参考资料与导航

实际资料核验与编译/运行记录位于 writing-plans/computer-graphics/evidence/research-31.md、31-cpu.txt。数值范围只对应当前质量、刚度、初态和时长。

系列入口 · 上一篇：骨骼怎样带动网格 · 下一篇：固定步长与简单碰撞。