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

Question

我正在使用不同的数值方法来解决阻尼简谐振子的常微分方程（ODE）。我通过计算预测能量与解析解之间的平均绝对误差（MAE）来比较每种积分方法的表现：
$$ MAE = \frac{1}{n} \sum_{i=0}^{n}\left |y_{analytical}-y_{numerical}\right| $$

然后，我对不同的时间步长计算得到的MAE，并将结果绘制成对数坐标系的图，如下所示。

log (MAE) vs. log(Time_step)

MAE和时间步长之间的关系符合我的预期（Verlet方法的误差随时间步长的平方增长，而Euler-Cromer方法的误差随时间步长线性增长），但我注意到Verlet方法在大约10^(-4)秒时有一个转折点。这似乎有点大，我原本期待转折点出现在更接近10^(-8)秒的时间步长，因为我使用的是numpy的float64类型，这样大约有15到17位小数的精度。

我接着绘制了每个时间步长的最大和最小误差（排除了第0次迭代，因为那是初始条件，对于数值方法和解析方法都是相同的），结果如下： log (Min Err) vs. log(Time_step) 和 log (Max Err) vs. log(Time_step)

再次绘制最大误差时，我得到的最佳时间步长与之前的图相似，但绘制最小误差时（这些误差总是在初始条件后的前几次迭代中出现），我发现误差在10^(-4)秒时似乎趋于平稳，接近10^(-15)焦耳的能量误差。

由于最小误差的平稳化，进一步增加时间步长到10^(-4)秒似乎不会提高Verlet方法的精度，但我无法解释为什么最大误差在这个点之后会增加。

我想到的一个解释是，由于float64的舍入误差，当数值达到大约10^(-16)时会出现这种情况。我手动检查了运行Verlet方法后得到的位置、速度和加速度，它们的最小值大约在10^(-9)的数量级，远离10^(-16)。

(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

数值方法常微分方程舍入误差阻尼简谐运动 verlet方法时间步长平均绝对误差能量误差

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

1 个回答

撰写回答