Verlet 方法求解阻尼简谐运动 ODE 的最佳时间步长

Question

我使用不同的数值方法来理解不同类型的积分器在不同时间步长产生的结果。我通过计算预测能量的平均绝对误差与解析解来比较每种积分方法的性能：

$$ MAE = \frac{1}{n} \sum_{i=0}^{n}\left |y_{analytical}-y_{numerical}\right| $$

然后，对于不同的时间步长，我计算所得的 MAE 并将结果绘制在对数与对数图中，如下所示。

MAE 和时间步之间的关系符合我的预期（Verlet 方法按二次方缩放，Euler Cromer 方法按线性缩放），但我注意到 Verlet 方法在大约 10^(-4) s 处有一个转折点。这看起来有点太大了，我希望在接近 10^(-8) s 的时间步长处出现转折点，因为我使用的是 numpy 的 float64，因此精度大约有 15 到 17 位小数。

我绘制了每个时间步长获得的最大和最小误差（不包括迭代 0，因为这些是初始条件，对于数值方法和分析方法来说都是相同的），结果如下： 和

再次绘制最大误差时，与之前的图相比，我获得了相似值的最小值，但绘制了获得的最小误差（这些总是发生在初始条件后的前几次迭代中），我发现误差似乎在以下位置变平： 10^(-4) s，接近误差约为 10^(-15) J 的能量。

由于最小误差的扁平化，进一步超过 10^(-4) s 并不会提高 Verlet 方法的精度是有道理的，但我无法解释为什么最大误差在这一点之后会增长。

我想到的一个解释是 float64 引起的舍入误差，当值达到大约 10^(-15)、10^(-17) 时，就会发生这种误差。我手动检查了运行 verlet 方法产生的位置、速度和加速度，但它们的最低值是 10^(-9)，与 10^(-15) 相差甚远。

(1) 当我根据分析方法和 Verlet 方法计算残差时，是否可能会引入舍入误差？

(2)还有其他更合适的误差计算方法吗？ （我认为 MAE 很合适，因为 verlet 方法会围绕真实系统值振荡）

(3) 是否可以进行一些调整来显示我的分析中可能存在的缺陷，我已经广泛地查看了我的代码，但我无法找到任何错误，此外，我编码的 Verlet 方法确实存在一个按二次方缩放的错误时间步长让我认为代码本身很好。 （也许可能的尝试是使用 float128 并确保在所有计算中使用它，然后看看上面的图是否不同？）

预先感谢您对上述问题的任何帮助

这是我计算 MAE 的方法：

def calculate_mean_absolute_error(numerical_method_energy,
analytical_method_energy): 
    residual = np.abs(analytical_method_energy - numerical_method_energy)

    mean_absolute_error = sum(residual) / len(residual)
    
    return mean_absolute_error

Answer 1

在比较一阶方法和二阶方法时，这完全符合预期。

计算误差有两个影响，实际计算的数值步长误差和浮点误差。每单位步长的理论误差为 ~h^p，而每单位步长的浮点误差与步数 1/h 成正比。

V 形的钩子或扭结是这些误差平衡的地方，因此在 h ~ mu^(1/(p+1)) 处。对于二阶方法，该值为 10^(-5)；对于使用标准双精度类型的一阶方法，该值为 10^(-8)。

更多详细信息请参见https://math.stackexchange.com/questions/3609842/numerical-analysis-heuns-method