C++实现欧拉法求解弹簧振子运动方程：数值发散过快问题

核心原因：显式欧拉法的稳定性阈值 #

对于质量m=1、弹簧系数k=1的弹簧振子，运动方程为$\ddot{x} = -x$。显式欧拉法求解这类二阶线性保守系统时，存在严格的稳定步长限制：只有当迭代步长$\Delta t < 2\sqrt{\frac{m}{k}}$（此处即$\Delta t < 2$）时，数值解才会呈现缓慢的能量增长；若步长等于或超过2，数值解会直接进入指数发散状态，迭代几次就会出现振幅爆炸。

代码层面的常见错误点 #

1. 步长设置超标：如果代码中$\Delta t$设为2或更大（比如3、5），会直接触发不稳定条件，10次迭代内必然严重发散。
2. 符号错误：若将加速度项写成$+x$而非$-x$，运动方程会从简谐振动变成指数增长的不稳定系统，数值解会瞬间发散。
3. 递推逻辑错误：显式欧拉法要求用当前时刻的位置$x_n$和速度$v_n$计算下一时刻的值，正确的递推式应为：

   double x_next = x_current + v_current * dt;
   double v_next = v_current - x_current * dt; // 因m=1,k=1，简化为-x_current*dt
   

   若错误地先用更新后的速度计算位置，或者颠倒了递推顺序，也可能加剧发散速度。

解决方向 #

• 立即将步长$\Delta t$调整到2以内（比如0.1、0.5），此时数值解会呈现预期的缓慢能量增长趋势；
• 改用半隐式欧拉、Verlet积分等对保守系统更稳定的数值方法，这类方法无严格步长稳定限制（或限制宽松），能长期保持能量近似守恒。

内容的提问来源于stack exchange，提问作者ugur