最近微分方程课上有一个关于微分方程数值解的小组作业,这里记录一下微分方程数值解的稳定性问题。关于一个微分方程 $x^\prime(t)=f(t,x)$,给定初值 $x(0)=x_0$,称为一个初值问题。这个作业使用的 Euler 方法的想法可以看作是 Taylor 展开

取了一阶近似然后累加的结果。Euler 法的意思就是求一个数组,让 $y_{k+1}=y_k + \Delta t f(t_k, y_k)$ 不断累加,就可以取到 $y(t)$ 的近似解。

然而由于 Euler 方法比较粗糙,得出的数值解就不一定收敛。尤其是我们研究的例子 $y^\prime=e^t \sin y$。这个方程有趣的地方是,它本身是可以给出解析解的,解就是 $\DeclareMathOperator{\arccot}{arccot}y=2\arccot(e^{-e^t+c})$,并且当 $t$ 趋于无穷的时候,$y$ 是收敛到 $\pi$ 的,但 Euler 方法在这个方程上却不收敛。

对于这个问题,A.B.H 大佬给出了一个漂亮的解释。现在我们关心的是 Euler 法计算出的计算值 $y_{k}$ 与理论值 $Y_k=y(t_k)$ 之间的绝对误差是怎么变化的。那么就设 $E_k=Y_k-y_k$。

而对于这个误差的估计,是有

这样的一个递归式。然后再一次对 $Y(t_k+\Delta t)$ 做一个 Taylor 展开。

然后运用一下中值定理,

就能有一个很好的对误差的估计,

终于就可以估计误差了!但是还要经过一个非常痛苦的计算可以算出来 $y^\prime(t)=\frac{2e^{c+e^t+t}}{e^2c+e^{2e^t}}$ 和 $y^{\prime\prime}(t)=\frac{2(e^t+1)e^{c+e^t+t}}{e^{2c}+e^{2e^t}}-\frac{4e^{c+3e^t+2t}}{(e^2c+e^{2e^t})^2}$。当然这个不得不让 Mathematica 代劳了。从这两个公式的分子和分母的函数的阶可以看出来当 $t \rightarrow \infty$ 的时候函数值都是趋于零的。

于是在忽略二阶导数的误差估计中

于是就可以看出当 $t\rightarrow\infty$ 的时候,就有一个 $e^{t_k}$ 这样一个非常大的误差。所以 Euler 方法最终会不收敛。

这个解法是我和 A.B.H 大佬在清华学堂二楼讨论了一个小时的结果,最后还是 A.B.H 给出了关键步骤的想法。大概 teamwork 被带飞就是这种体验吧。最后是一张盗图,当时的讨论用的黑板的凌乱的公式。