如何在p5.js中模拟波动方程，解决振幅异常增大问题

错误点与修正方案 #

• 离散公式逻辑错误
  你给出的离散形式代码中对波动方程的转换是错误的：波动方程右侧(d²u/dx²)是位移对空间坐标x的二阶偏导，离散后的分母应该是空间采样步长的平方dx²，你错误地使用了时间步长作为分母。即使你前面的积分代码中除以了d2x，也需要先确认d2x的取值确实等于两个相邻粒子x方向间距的平方，否则空间二阶导的计算从根本上就是错的。

• 不满足CFL稳定性条件
  显式时间积分求解波动方程有硬性的稳定约束，也就是CFL条件：c * dt / dx ≤ 1，其中c是你设定的波速，dt是时间步长，dx是相邻粒子的x方向间距。只要这个条件不满足，数值解必然会出现振幅无限制放大的不稳定问题，和你当前的现象完全吻合。你可以调小波速、调小时间步长，或者调大粒子的x方向间距来满足这个条件。

• 积分格式选择不当
  你使用的前向欧拉积分格式天生不适合二阶波动方程的模拟，数值稳定性差，建议换成专门适配波动场景的蛙跳积分格式：

// 蛙跳积分示例
vel.y += acc.y * dt / 2; // 先更新半时间步的速度
pos.y += vel.y * dt;
// 下一帧计算新的acc之后，再更新另一半速度
vel.y += newAcc.y * dt / 2;

• 缺少阻尼抑制数值误差
  就算满足了所有稳定条件，数值计算的截断误差也可能慢慢累积导致振幅缓慢上升，可以加一个极小的阻尼项抑制异常增长，不会对波的正常传播产生明显影响：

vel.y *= 0.9995; // 每次速度更新后乘以接近1的阻尼系数

内容的提问来源于stack exchange，提问作者Álvaro Rodrigo