常微分方程的数值解

常微分方程的离散化

一阶微分方程的一般形式为:

\[\begin{cases} \frac{dy}{dx} = f(x,y)&a\leq x\leq b \\ y(a) = y_0 \tag{1} \end{cases} \]

假定\(f(x,y)\)连续，且关于\(y\)满足\(Lipschitz\)条件，即存在常数L，使得\(|f(x,y)-f(x,\overline{y})|\leq L|y-\overline{y}|\)
由常微分方程理论知，初值问题(1)的解存在且唯一。
所谓数值解法，就是求问题(1)的解\(y(x)\)在若干点

\[a = x_0 < x_1 < x2 < ...<x_N = b \]

处的近似值\(y_n(n=1,2,...,N)\)的方法。
\(y_n(n=1,2,...,N)\)称为问题(1)数值解,\(h_n=x_{n+1}-x_n\)的步长。

建立数值解法，首先要将微分方程离散化，采用如下三种方法：

（1）差商近似导数

用向前差商\(\frac{y(x_{n+1})-y(x_n)}{h}\)代替\(y^{'} (x_n)\)带入微分方程(1)得到：

\[\frac{y(x_{n+1})-y(x_n)}{h} \approx f(x_n,y(x_n)) \]

化简得$$y(x_{n+1}) \approx y(x_n)+hf(x_n,y(x_n))$$
用\(y(x_n)\)的近似值\(y_n\)代入上式右端，所得结果作为\(y(x_{n+1})\)的近似值，记为\(y_{n+1}\),得到：

\[y_n+1 = y_n + hf(x_n,y_n) \quad (n=0,1,...)\tag{2} \]

这样，问题(1)的近似解可通过求解以下问题得到

\[\begin{cases} y_{n+1} = y_n+hf(x_n,y_n) & (n=0,1,...)\\ y_0 = y(a) \tag{3} \end{cases} \]

用步长和前面一个点的斜率计算后面一个点。

（2）数值积分方法

将问题(1)表示成积分形式，用数值积分方法离散化，如：

\[y(x_{n+1})-y(x_n) = \int^{x_{n+1}}_{x_n}f(x,y(x))dx \quad (n=0,1,...) \tag{4} \]

右边积分用矩形公式或者梯形公式近似计算。

Euler方法以及改进的Euler方法

1.Euler方法

由公式(3):

\[\begin{cases} y_{n+1} = y_n+hf(x_n,y_n) & (n=0,1,...)\\ y_0 = y(a) \end{cases} \]

算出\(y(x_n)\)的近似值\(y_n\),这组公式称为向前Euler公式

如果用向后差商代替导数，即

\[y^{'}(x_{n+1})\approx \frac{y(x_{n+1})-y(x_n)}{h} \]

则得到：

\[\begin{cases} y_{n+1} = y_n+hf(x_{n+1},y_{n+1}) & (n=0,1,...)\\ y_0 = y(a) \tag{5} \end{cases} \]

向后Euler公式属于隐式公式，实际计算要复杂得多，一般用迭代法求解，迭代公式为：

\[\begin{cases} y_{n+1}^{0} = y_n+hf(x_{n},y_{n})\\ y_{n+1}^{k+1} = y_n+hf(x_{n+1},y_{n+1}^{k}) & (k=0,1,2,...) \end{cases} \tag{6} \]

迭代过程如下：

用显式公式求出\(y_{n+1}\)作为\(y^{(0)}_{n+1}\),即

\[y^{(0)}_{n+1} = y_n + hf(x_n,y_n)\\ y^{(1)}_{n+1} = y_n + hf(x_{n+1},y_{n+1}^{0})\\ y^{(2)}_{n+1} = y_n + hf(x_{n+1},y_{n+1}^0)\\ \vdots\\ y_{n+1}^{k+1} = y_n+hf(x_{n+1},y_{n+1}^{k}) \]

当\(|y_{n+1}^{k+1}-y_{n+1}^{k} < \varepsilon|\)时，取\(y_{n+1} = y^{k+1}_{n+1}\)

Euler方法得误差估计

对于向前Euler公式，右端的\(y_n\)都是近似的，用它计算\(y_{n+1}\)会有累积误差，分析累积误差比较复杂，先讨论比较简单的局部截断误差。
假定（3）右端的\(y_n\)没有误差，即\(y_n=y(x_n)\),由此得到

\[y_{n+1} = y_n+hf(x_{n},y_{n})\tag{7} \]

由Taylor展开得到精确值\(y(x_{n+1})\)是

\[y(x_{n+1}) = y(x_n) + hy^{'}(x_n)+\frac{h^2}{2}y^{''}(x_n)+O(h^3)\tag{8} \]

(7)、(8)相减得到：

\[y(x_{n+1})-y_{n+1} = \frac{h^2}{2}y^{''}(x_n)+O(h^3) \]

即局部截断误差是\(h^2\)阶的
数值算法精度定义为：

若一种算法截断误差为\(O(h^{p+1})\),则称该算法具有\(p\)阶精度。
显然p越大精度越高，向前Euler方法是一阶的，因此精度不高。

改进的Euler方法

（1）梯形公式
利用数值积分方法将微分方程离散化，若用梯形公式（4）计算右端积分，即

\[\int^{x_{n+1}}_{x_n}f(x,y(x))dx\approx\frac{h}{2}[f(x_n,y(x_n))+f(x_{n+1},y(x_{n+1}))] \]

依然用近似值代替，得到计算公式

\[y_{n+1}=y_n+\frac{h}{2}[f(x_n,y(x_n))+f(x_{n+1},y(x_{n+1}))] \]

这是求初值问题(1)的梯形公式，向前Euler公式为矩形公式。梯形公式为二阶方法。

梯形公式也是隐式格式，需要迭代求解，迭代公式

\[\begin{cases} y_{n+1}^{0} = y_n+hf(x_{n},y_{n})\\ y_{n+1}^{k+1} = y_n+\frac{h}{2}[f(x_n,y_n)+f(x_{n+1},y_{n+1}^{k})] & (k=0,1,2,...) \end{cases} \]

函数\(f(x,y)\)满足\(Lipschitz\)条件，容易看出

\[|y_{n+1}^{(k+1)}-y_{n+1}^{(k)}|\leq\frac{hL}{2}|y_{n+1}^{(k)}-y_{n+1}^{k-1}| \]

因此，当\(0<\frac{hL}{2}<1\)时，迭代收敛。
但是这样做计算量大，因此导出一种新的方法————改进的Euler方法。

(2)改进的Euler方法

将Euler公式与梯形公式结合使用：先用Euler公式求出\(y_{n+1}\)的一个初步近似值\(\overline{y}_{n+1}\),称为预测值，然后用梯形公式校正求得近似值\(y_{n+1}\)，
即：

\[\begin{cases} \overline{y}_{n+1} = y_n+hf(x_{n},y_{n}) \quad预测\\ y_{n+1} = y_n+\frac{h}{2}[f(x_n,y_n)+f(x_{n+1},\overline{y}_{n+1})] \quad校正 \end{cases} \tag{6} \]

该预测--校正系统也叫改进的Euler方法。
该方法为二阶方法。

为便于上机编写程序，式（11）常改写成

\[\begin{cases} y_p = y_n+hf(x_n,y_n)\\ y_q = y_n + hf(x_n+h,y_p)\\ y_{n+1}=\dfrac{1}{2}(y_p+y_q) \end{cases} \tag{12} \]

龙格-库塔（Runge-Kutta）方法

Eular方法的基本思想是用差商代替导数，由微分中值定理，有

\[\dfrac{y(x_{n+1})-y(x_{n})}{h} = y^{'}(x_n+ \theta h),\quad 0<\theta < 1 \]

由方程\(y^{'}=f(x,y)\),有

\[y(x_{n+1}) = y(x_n) + hf(x_n+\theta h,y(x_n+\theta h)) \tag{13} \]

记\(\overline{K} = f(x_n+\theta h,y(x_n+\theta h))\),
称为区间\([x_n,x_{n+1}]\)上的平均斜率。
给出一种斜率\(\overline{K}\),就可根据上式导出一种算法。

向前Euler公式，取\(\overline{K} = f(x_n,y_n)\),精度很低。

改进Euler公式，取\(\overline{K} = \dfrac{1}{2}[f(x_n,y_n)+f(x_{n+1},\overline{y}_{n+1})]\),提高了精度。

如果我们在区间\([x_n,x_{n+1}]\)内多取几个点，将其斜率加权平均作为\(\overline{K}\)，自然可以构造出精度更高的计算公式，这就是龙格-库塔方法的基本思想。

二阶龙格-库塔公式

在区间\([x_n,x_{n+1}]\)内，除\(x_n\)外再取区间中任一点\(p\),则有：

\[x_{n+p}=x_n+ph \quad 0 < p \leq 1 \]

取\(x_n\)和\(x_{n+p}\)两点的斜率\(k_1\)和\(k_2\)加权平均作为平均斜率，可构造出以下计算格式：

\[y_{n+1} = y_n+h[(1-\lambda)k_1+ \lambda k_2] \]

其中:

\[k_1 = f(x_n,y_n)\\ k_2 = f(x_{n+p},y_n+phk_1) \]

该格式为二阶的Runge-Kutta格式，其中\(p\),\(\lambda\)为待定参数。在确定\(p\),\(\lambda\) 的同时，使二阶Runge-Kutta格式具有较高的精度，具体推导如下：

假设\(y_n = y(x_n)\)，有：

\[k_1 = f(x_n,y_n) = y^{'}(x_n)\\ k_2 = f(x_{n+p},y_n+phk_1) = f(x_n+ph,y_n+phk_1) \]

将\(k_2\)在点\((x_n,y_n)\)处展开为Taylor级数：

\[\begin{aligned} k_2 &= f(x_n,y_n)+phf_x(x_n,y_n)+phy^{'}(x_n)f_y(x_n,y_n)+O(h^2)\\ &= y^{'}(x_n) + phy^{''}(x_n)+O(h^2) \end{aligned} \]

将\(k_1\)和\(k_2\)代入\(y_{n+1}=y_n+h[(1-\lambda)k_1+\lambda k_2]\)
得到：

\[y_{n+1} = y(x_n)+hy^{'}(x_n)+ \lambda ph^2y^{''}(x_n)+O(h^3) \]

再将\(y(x_{n+1})\)在点\((x_n,y_n)\)处展开为Taylor级数：

\[y(x_{n+1}) = y(x_n+h) = y(x_n)+h^2y{''}(x_n)+O(h^3) \]

比较两式，得：

\[y(x_{n+1})-y_{n+1} = (\dfrac{1}{2} h^2-\lambda ph^2)y{''}(x_n)-\lambda O(h^3) \]

因此，当\(\lambda p = \dfrac{1}{2}\)时，二阶Runge-Kutta格式具有二阶精度。二阶Rungfe-Kutta格式有很多具体形式：

若取\(\lambda = \frac{1}{2},p = 1\)，则可得到改进的Euler格式；
若取\(\lambda = \frac{3}{4}, p = \frac{2}{3}\),则可得Heun(休恩)格式；
若取\(\lambda = 1, p = \frac{1}{2}\), 则可得变形的Euler格式，即中点格式：

\[\begin{cases} y_{n+1} = y_n+hk_2\\ k_1 = f(x_n,y_n)\\ k_2=f(x_{n+\frac{1}{2}},y_n+\frac{h}{2}k_1) \end{cases} \]

四阶龙格-库塔公式

在区间\([x_n,x_{n+1}]\)内，使用两个不同的点可以构造出二阶Runge-Kutta格式。依此规律，在区间\([x_n,x_{n+1}]\)内，取3个不同的点可以构造出三阶Runge-Kutta格式；取4个不同的点可以构造出四阶Runge-Kutta格式。在实际中，应用最广泛的是四阶经典的Runge-Kutta格式：

\[\begin{cases} y_{n+1} = y_n + \cfrac{h}{6}(k_1 + 2k_2 + 2k_3+k_4)\\ k_1 = f(x_n, y_n)\\ k_2 = f(x_n+\frac{h}{2},y_n+\frac{h}{2}k_1)\\ k_3 = f(x_{n+\frac{1}{2}},y_n + \frac{h}{2}k_2)\\ k_4 = f(x_n + h,y_n+hk_3) \end{cases} \]

该方法是一种高精度的单步法，可达到四阶精度\(O(h^5)\)。

四阶Runge-Kutta法除经典的格式外，还有其它格式。比如Kutta格式和Gill格式。

线性多步法

以上介绍的数值解法都是单步法，因为它们在计算\(y_{n+1}\)时，都只用到前一步\(y_n\)的值，单步法的一般形式是：

\[y_{n+1} = y_n + h\phi(x_n,y_n,h) \quad n = 0,1,\dots,N-1 \]

其中，\(phi(x_n,y_n,h)\)称为增量函数。如何通过较多地利用前面的已知信息，如\(y_n,y_{n-1},\dots\)来构造高精度的算法，计算\(y_{n+1}\),这就是多步法的基本思想。

经常使用的是线性多步法。设\(r = 1\),即用\(y_n,y_{n-1}\)计算\(y_{n+1}\)

首先用数值积分方法离散化方程

\[y(x_{n+1}) - y(x_n) = \int _{x_n} ^ {x_{n+1}} f(x,y(x))dx \]

记\(f(x_n,y_n) = f_n\),\(f(x_{n-1},y_{n-1}) = f_{n-1}\),式中被积函数
\(f(x,y(x))\)用二节点\((x_{n-1},f_{n-1}), (x_n,f_n)\)的插值公式得到。
\(x >= x_n\)所以是外插。

\[\begin{aligned} f(x,y(x)) &= f_n\frac{x-x_{n}}{x_n-x_{n-1}} + f_{n-1}\frac{x-x_n}{x_{n-1}-x_n}\\ &= \frac{1}{h}[(x-x_{n-1})f_n-(x-x_n)f_{n-1}] \end{aligned} \]

在区间\([x_n,x_{n+1}]\)上积分可得

\[\int _{x_n} ^{x_{n+1}} f(x,y(x))dx = \frac{3h}{2}f_n - \frac{h}{2}f_{n-1} \]

于是\(y_{n+1} = y_n + \dfrac{h}{2}(3f_n-f_{n-1}) \qquad (20)\)
公式(20)截断误差为\(O(h^3)\)，比向前Euler公式精度提高一阶。

类似可推导\(r = 2,3,\dots\)时的公式，数值越大精度越高。

一阶微分方程组与高阶微分方程的数值解法

一阶微分方程组的数值解法

设有一阶微分方程组的初值问题

\[\begin{cases} y^{'} = f_i(x,y_1,y_2,\dots,y_m)\\ y_i(a) = y_{i0} \end{cases} (i = 1,2,\dots,m) \tag{21} \]

若记 \(y = (y1, y2,\dots,y_m )^T,y_0 = (y_{10},y_{20},\dots,y_{m_0})^{T},f = (f_1,f_2,\dots,f_m)^T\),则初值问题(21)可以写成向量的形式

\[\begin{cases} y^{'} = f(x,y)\\ y(a) = y_0 \end{cases} \tag{22} \]

如果向量值函数\(f(x,y)\)在区域D:\(a\leq x \leq b, y\in R^m\)连续，且关于\(y\)满足\(Lipschitz\)条件，即存在\(L > 0\),使得对\(\forall x \in [a,b],y_1,y_2 \in R^m\),都有

\[\left || da \right || \leq L \left || y_1 - y_2 \right|| \]

那么问题(22)在[a,b]上存在唯一解 $ y = y(x)$
问腿(22)形式上与(1)完全相同，故(1)的各种数值方法也适用于(22)

高阶微分方程的数值解法

简单来说， 高阶微分方程的初值问题可通过变量代换转化为一阶微分方程组的初值问题。
设有m阶常微分方程初值问题:

\[\begin{cases} y^(m) &= f(x,y,y^{'},\dots,y^{(m-1)}) & a \leq x \leq b \\ y(a) & = y_0, y^{'}(a) = y_0^{(1)},\dots,y^{(m-1)}(a) = y_0^{(m-1)} \end{cases} \]

引入新变量 \(y_1 = y, y_2 = y^{'},\dots,y_m = y^{(m-1)}\),问题(23)就化为一阶微分方程的初值问题。

刚性问题

在化学工程及自动控制等领域中，所涉及的常微分方程组初值问题常常是所谓的“刚性”问题。具体地说，对一阶线性微分方程组

\[\dfrac{dy}{dx} = Ay + \phi(x) \]

式中：\(y,\phi \in R^m\);A为m阶方阵。若矩阵A的特征值\(\lambda_i(i=1,2,\dots,m)\)满足关系

\[Re \lambda_{i} < 0, i = 1,2,...,m\\ \mathop{max} \limits _{1\leq i \leq m}|Re \lambda_i| \gg \mathop{min} \limits _{1\leq i \leq m}|Re \lambda_i| \]

则称方程组为刚性方程组或Stiff方程组。
称数

\[s = \mathop{max} \limits _{1\leq i \leq m}|Re \lambda_i|/ \mathop{min} \limits _{1\leq i \leq m}|Re \lambda_i| \]

为刚性比。理论上的分析表明，求解刚性问题所选用的数值方法最好是对步长h不做任何限制。