常微分方程的数值解

导语

本文介绍了常微分方程的初值问题、边值问题与本征值问题的数值解法，正确理解这些方法，对MATLAB或其它编程语言计算常微分方程有着极大帮助，同时而言也有利于后续偏微分方程的学习。

常微分方程回顾

首先再来认识一下微分方程和常微分方程。

微分方程是未知函数及其导数或微分的关系，例如：$\frac{\mathrm{d}y}{\mathrm{d}x}=2x$ 等。函数方程是关于未知数$x$的，而微分方程则是关于$x$的导数或微分的，这就是微分方程 (Differential Equation)。

常微分方程则是一个未知函数与其导数或微分的关系，“常”(Ordinary)便表示“一个”未知函数这种平常的形式 (平常是相较于多个未知函数的偏微分而言)。因此，这样的只存在一个未知函数的导数/微分的关系式称作常微分方程(Ordinary Differential Equation -> ODE)。

一般形式$n$阶常微分方程为 $F(x,y,y^{\prime },\dots ,y^{(n)})=0$，一阶常微分方程的形式为 $y^{\prime }(x)=\frac{\mathrm{d}y}{\mathrm{d}x}=k(x,y)$。

常微分方程的两大类问题为初值问题和边值问题。初值问题：$y$在$t=0$时的值是给定的，求解 $t>0$ 时方程的解；边值问题：方程在某些点的值被给定，求解其他点的值。

为了方便讨论，我们以一阶常微分方程为例，并加入定解条件 $y(x_0)=y_0$，因此得到方程组 $\{\begin{array}{l}{y}^{\prime }(x)=\frac{dy}{dx}=k(x,y)\\ y\left(x_0\right)=y_0\end{array}$ ，那么有了这样的方程组以后，如何求解$y(x)$呢？

微分方程有两种常见的物理图像，一种为对时间的微分，例如位移的导数为速度 $x^{\prime }=\frac{\mathrm{d}x}{\mathrm{d}t}$等；而另一种即为曲线的斜率 $y^{\prime }=\frac{\mathrm{d}y}{\mathrm{d}x}$，在此处使用了第二种物理图像进行推导。
$k(x,y)$表示平面上各点处的斜率。

欧拉法

上述方程组给出了平面上点的坐标以及点对应的斜率值，可以使用欧拉法对 $y(x)$ 进行迭代的近似求解。求解分析如下：
已知 $x=x_0$ 时 $y=y_0$ ，且斜率为 $y^{\prime }\left(x_0\right)=k\left(x_0,y_0\right)$ ，那么就可以用过该点且斜率为 $y^{\prime }\left(x_0\right)$ 的直线近似求解 $x_0$ 附近一点 $x_1$ 的函数值。

构建直线：$y(x)=y_0+k\left(x_0,y_0\right)\left(x-x_0\right)$，那么，在 $x=x_1$ 时，求得
$$y\left(x_1\right)=y_0+k\left(x_0,y_0\right)\left(x_1-x_0\right)$$
即在 $x=x_1$ 时，$y=y\left(x_1\right)$，且斜率为 $y^{\prime }\left(x_1\right)=k\left(x_1,y_1\right)$ 。重复以上步骤，可以求得 $y\left(x_2\right),y\left(x_3\right),\cdots ,y\left(x_n\right)$

如果等间隔取 $x$ 值，即 $x_1-x_0=x_2-x_1=x_3-x_2=\cdots =x_n-x_{n-1}=h$ ，也即 $x_n=x_0+n\cdot h$ 可以得到
$$\begin{array}{c}y\left(x_1\right)=y_0+k\left(x_0,y_0\right)\cdot h\\ y\left(x_2\right)=y_1+k\left(x_1,y_1\right)\cdot h\\ ⋮\\ y\left(x_n\right)=y_{n-1}+k\left(x_{n-1},y_{n-1}\right)\cdot h\end{array}$$
总结为：
$$\{\begin{array}{l}y\left(x_n\right)=y_{n-1}+k\left(x_{n-1},y_{n-1}\right)\cdot h\\ x_n=x_0+n\cdot h\end{array}$$
该公式便是欧拉公式。

推导结束之后我们再来找找欧拉法的bug。可以注意到，欧拉法最大的问题在于斜率的选取，公式 $y\left(x_n\right)=y_{n-1}+f\left(x_{n-1},y_{n-1}\right)\cdot h=y_{n-1}+y^{\prime }(x_{n-1})\cdot h$ 中，我们用$x_{n-1}$处的斜率$y^{\prime }(x_{n-1})$乘间距$h$表示 $y(x_n)$ 与 $y(x_{n-1})$ 的差值，这样显然存在不小的误差。

如果读者是从按专栏顺序读下来的，那应该可以体会到，数值计算就是妥协的艺术，欧拉法便是在此处对斜率做了比较大的妥协。而我们可以使用更好的“妥协方法”将斜率表示出来，这样的妥协艺术就是 后退欧拉法，梯形法，改进欧拉法。

欧拉法：$y_{n+1}=y_n+hk\left(x_n,y_n\right)$；
后退欧拉法：$y_{n+1}=y_n+hk\left(x_{n+1},y_{n+1}\right)$；
梯形法：$y_{n+1}=y_n+\frac{h}{2}\left(k\left(x_n,y_n\right)+k\left(x_{n+1},y_{n+1}\right)\right)$；
改进欧拉法：$y_{n+1}=y_n+\frac{h}{2}$