常微分方程的数值解

导语

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

常微分方程回顾

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

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

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

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

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

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

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

欧拉法

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

构建直线： y ( x ) = y 0 + k ( x 0 , y 0 ) ( x − x 0 ) y(x)=y_{0}+k\left(x_{0}, y_{0}\right)\left(x-x_{0}\right) y(x)=y0​+k(x0​,y0​)(x−x0​)，那么，在 x = x 1 x=x_{1} x=x1​ 时，求得 y ( x 1 ) = y 0 + k ( x 0 , y 0 ) ( x 1 − x 0 ) y\left(x_{1}\right)=y_{0}+k\left(x_{0}, y_{0}\right)\left(x_{1}-x_{0}\right) y(x1​)=y0​+k(x0​,y0​)(x1​−x0​) 即在 x = x 1 x=x_{1} x=x1​ 时， y = y ( x 1 ) y=y\left(x_{1}\right) y=y(x1​)，且斜率为 y ′ ( x 1 ) = k ( x 1 , y 1 ) y^{\prime}\left(x_{1}\right)=k\left(x_{1}, y_{1}\right) y′(x1​)=k(x1​,y1​) 。重复以上步骤，可以求得 y ( x 2 ) , y ( x 3 ) , ⋯   , y ( x n ) y\left(x_{2}\right), y\left(x_{3}\right), \cdots, y\left(x_{n}\right) y(x2​),y(x3​),⋯,y(xn​)

如果等间隔取 x x x 值，即 x 1 − x 0 = x 2 − x 1 = x 3 − x 2 = ⋯ = x n − x n − 1 = h x_{1}-x_{0}=x_{2}-x_{1}=x_{3}-x_{2}=\cdots=x_{n}-x_{n-1}=h x1​−x0​=x2​−x1​=x3​−x2​=⋯=xn​−xn−1​=h ，也即 x n = x 0 + n ⋅ h x_{n}=x_{0}+n \cdot h xn​=x0​+n⋅h 可以得到 y ( x 1 ) = y 0 + k ( x 0 , y 0 ) ⋅ h y ( x 2 ) = y 1 + k ( x 1 , y 1 ) ⋅ h ⋮ y ( x n ) = y n − 1 + k ( x n − 1 , y n − 1 ) ⋅ h \begin{gathered} 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 \\ \vdots \\ y\left(x_{n}\right)=y_{n-1}+k\left(x_{n-1}, y_{n-1}\right) \cdot h \end{gathered} y(x1​)=y0​+k(x0​,y0​)⋅hy(x2​)=y1​+k(x1​,y1​)⋅h⋮y(xn​)=yn−1​+k(xn−1​,yn−1​)⋅h​ 总结为： { y ( x n ) = y n − 1 + k ( x n − 1 , y n − 1 ) ⋅ h x n = x 0 + n ⋅ h \left\{\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}\right. {y(xn​)=yn−1​+k(xn−1​,yn−1​)⋅hxn​=x0​+n⋅h​ 该公式便是欧拉公式。

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

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

欧拉法： y n + 1 = y n + h k ( x n , y n ) y_{n+1}=y_{n}+h k\left(x_{n}, y_{n}\right) yn+1​=yn​+hk(xn​,yn​)； 后退欧拉法： y n + 1 = y n + h k ( x n + 1 , y n + 1 ) y_{n+1}=y_{n}+h k\left(x_{n+1}, y_{n+1}\right) yn+1​=yn​+hk(xn+1​,yn+1​)； 梯形法： y n + 1 = y n + h 2 ( k ( x n , y n ) + k ( x n + 1 , y n + 1 ) ) 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) yn+1​=yn​+2h​(k(xn​,yn​)+k(xn+1​,yn+1​))； 改进欧拉法： y n + 1 = y n + h 2 ( k ( x n , y n ) + k ( x n + 1 , y n + h k ( x n , y n ) ) ) y_{n+1}=y_{n}+\frac{h}{2}\left(k\left(x_{n}, y_{n}\right)+k\left(x_{n+1}, y_{n}+h k\left(x_{n}, y_{n}\right)\right)\right) yn+1​=yn​+2h​(k(xn​,yn​)+k(xn+1​,yn​+hk(xn​,yn​)))

龙格库塔法

龙格库塔法核心思路与欧拉法相似，甚至欧拉法就是龙格库塔法的一种

由中值定理，得 y ( x n + 1 ) − y ( x n ) h = y ˙ ( x n + θ k )   ( 0 < θ < 1 ) \frac{y(x_{n+1})-y(x_n)}{h}=\dot{y}(x_n+\theta k)\ (0