数值确定我计算的一阶 ODE 解是否稳定/不稳定的算法?

计算科学 数字 显式方法
2021-12-14 22:33:05

我们被分配了一项任务,我们必须确定 Dahlquist 方程的数值解x˙=λx, (λ=7) 时间步长0.5,0.25,0.125使用显式欧拉方法。当我对它进行编码并在 MATLAB 上绘制解决方案时,我看到解决方案在时间步长上非物理地振荡0.50.25. 我在考虑是否可以设计一种算法来预测解决方案对于任何一般 ODE 是否稳定x˙=f(x).

任何可以适用于上述类型的一般 ODE 的稳健算法的想法?

PS:我在想我也许可以看到数值解的斜率和解点处的解析导数是否具有相同的符号(即,如果xn+1xnδtf(xn)<0在任何时候(tn,xn) 意味着存在非物理振荡)但是当数值解和解析解都单调增加/减少时它不起作用(数值方法会随着时间的增加而变为非常大的正值或负值t)。

1个回答

首先,您可以为 ODE 步进器生成的数据定义可微插值,然后绘制残差

Δ(t):=y~(t)f(y~,t)

其中是解决方案骨架的插值。(在 matlab 中,不是解决方案骨架 pchip 的默认插值器吗?应该这样做;尽管三次 Hermite 插值可能会更好,因为很容易计算。)这几乎肯定会向您展示您的步进器没有做你想做的。如果您希望以编程方式确定是否使用不同的方法进行解析,我可能会计算y~y~(t)Δ

在 IMO 中,使用残差是您获取所需信息的最简单方法。但是,严格来说,我没有按照所述回答您的问题,因为您想知道是否有后验方法来确定步进器与 rhs的交互是否导致它不稳定。事实证明,如果步进器在初始值问题上不稳定,那么它在终值问题上通常是稳定的。您从 IVP 解决方案的最终位置及时向后求解方程。如果这两种解决方案不大致相同,那么您就知道步进器在一个或另一个方向上不稳定。我会使用残差来确定哪个是哪个。f