求解阻尼简谐运动ODE的Verlet方法最优时间步长问题咨询

我正在用不同数值方法求解阻尼简谐振动的常微分方程（ODE），通过计算预测能量与解析解的**平均绝对误差（MAE）**对比各积分方法的性能，MAE计算公式为：
$$ MAE = \frac{1}{n} \sum_{i=0}^{n}\left |y_{analytical}-y_{numerical}\right| $$

针对不同时间步长，我计算了对应的MAE并绘制了双对数（log vs. log）曲线（图：log(MAE)与log(时间步长)的关系）。MAE与时间步长的关系符合预期——Verlet方法呈二次缩放，Euler-Cromer方法呈线性缩放，但Verlet方法的MAE曲线在约$10^{{-4} \text{s}$处出现转折点，这比我的预期早得多：由于我使用numpy的float64精度（约15-17位十进制有效数字），原本预期转折点会在$10}{-8}\ \text{s}$附近。

随后我又绘制了各时间步长下的最大、最小误差曲线（排除初始条件对应的第0次迭代）：

• log(最小误差)与log(时间步长)的关系：最小误差均出现在初始条件后的前几次迭代，在$10^{{-4} \text{s}$处趋于平稳，能量误差接近$10}{-15}\ \text{J}$
• log(最大误差)与log(时间步长)的关系：最优时间步长的趋势和之前的MAE曲线相近，但时间步长小于$10^{-4}\ \text{s}$后，最大误差反而增大

最小误差的平稳性说明，当时间步长小于$10^{{-4} \text{s}$时，Verlet方法的精度无法再提升，但我无法解释为何此时最大误差会上升。我猜测可能是float64的舍入误差（通常在数值量级达到$10}{-16}$时出现），但手动检查Verlet方法得到的位置、速度、加速度值，其最低量级为$10^{{-9}$，远未达到$10}{-16}$的阈值。

疑问 #

• 计算解析解与Verlet方法结果的残差时，是否可能引入舍入误差？
• 是否存在更合适的误差计算方式？（我原本认为MAE适合Verlet方法，因其结果围绕真实值振荡）
• 我的分析是否存在可改进之处？我已反复检查代码未发现Bug，且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

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