数值分析学习笔记--常微分方程初值问题的数值解法

应用纳维-斯托克斯方程对导管中的气流进行模拟

如有疏漏，敬请指正。

对于一个常微分方程(ODE)

$$\begin{cases} \frac{\mathbf{d}y}{\mathbf{d}x}=f(x,y) \\ y(a)=y_0 \\ \end{cases} ,a \leq x \leq b$$

只有一些特殊形式的$f(x,y)$，才可以求它的解析解。对于大多数常微分方程的初值问题，只可以计算其数值解。

对区间$[a,b]$进行一个剖分，确定点列$\lbrace x_n,n=1,2, \cdots ,m \rbrace$，使用函数在这一组点列上的函数值$\lbrace y(x_n), n=1,2, \cdots ,m \rbrace$描述$y(x)$，求解$y(x_n)$的数值方法，是为有限差分方法。

约定$y(x_n)$是常微分方程准确解的值，$y_n$是$y(x_n)$的近似值。把区间分割为小区间的形式，称为数值离散方法。

除了有限差分方法外，有限元方法与有限体积方法也可以求解ODE的数值解。

欧拉公式

欧拉方法的图示，待求的曲线为蓝色，它的折线近似为红色。

对区间$[a,b]$进行等距剖分，记步长$h=\frac{b-a}{m}$，那么$\lbrace x_n=a+nh,n=0,1, \cdots ,m \rbrace$。可以使用差商近似导数来求解ODE。

• 向前差商
  作出$y’(x)$在$x=x_n$处的一阶向前差商 $$y'(x_n) \approx \frac{y(x_{n+1})-y(x_n)}{h}$$ 而由微分方程可得 $$y'(x_n)=f(x_n,y(x_n))$$ 这样就可以得到计算$y_{n+1}$的向前欧拉公式 $$y_{n+1}=y_n+hf(x_n,y_n)$$ 不难发现，向前欧拉公式是一个显式求解格式。
• 向后差商
  作出$y’(x)$在$x=x_{n+1}$处的一阶向后差商 $$y'(x_{n+1}) \approx \frac{y(x_{n+1})-y(x_n)}{h}$$ 又因为$y’(x_{n+1})=f(x_{n+1},y(x_{n+1}))$，得到向后欧拉公式 $$y_{n+1}=y_n+hf(x_{n+1},y_{n+1})$$ 不难发现，这是一个隐式求解格式。设$f(x,y)$对$y$满足Lipschitz条件，那么可以使用Picard格式求解向后欧拉公式。在实务中，一般先用显式格式先算出初始值，再用隐式公式进行进一步的迭代。这样的过程被称为预估-校正过程。向后欧拉公式的预估-校正式如下所示： $$\begin{cases} \bar{y}_{n+1}=y_n+hf(x_n,y_n) \\ y_{n+1}=y_n+hf(x_{n+1},\bar{y}_{n+1}) \\ \end{cases}$$
• 中心差商
  类似前面的方法，可得中心差商式为 $$y_{n+1}=y_{n-1}+2hf(x_n,y_n)$$ 这不是一个稳定的格式，因此在实务中不予采用。

欧拉公式的局部截断误差

对$y(x_{n+1})$在$x_n$进行泰勒展开

$$y(x_{n+1})=y(x_n)+hf(x_n,y(x_n))+\frac{h^2}{2!}y''(\zeta_n),x_n \leq \zeta_n \leq x_{n+1}$$

不难发现局部截断误差是

$$T_{n+1}=\frac{h^2}{2!}y''(\zeta_n)$$

欧拉公式的改进

对于方程$\frac{\mathbf{d}y}{\mathbf{d}x}=f(x,y)$在区间$[x_n,x_{n+1}]$上进行积分再移项

$$y(x_{n+1})=y(x_n)+\int_{x_n}^{x_{n+1}}f(x,y(x))\mathrm{d}x$$

取梯形公式计算数值积分

$$\int_{x_n}^{x_{n+1}}f(x,y) \mathrm{d}x=\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_n)+f(x_{n+1},y_{n+1}))$$

使用显式欧拉公式与梯形公式组成预估-校正公式

$$\begin{cases} \bar{y}_{n+1}=y_n+hf(x_n,y_n) \\ y_{n+1}=y_n+\frac{h}{2}(f(x_n,y_n)+f(x_{n+1},\bar{y}_{n+1})) \end{cases}$$

这被称为改进的欧拉公式。

Runge-Kutta方法

在二阶Runge-Kutta方法中，使用$f(x,y)$在点$(x_n,y(x_n))$以及其附近的点$(x_n+ah,y(x_n)+bhf(x_n,y(x_n)))$上的值的线性组合得到数值公式

$$y_{n+1}=y_n+h(c_1f(x_n,y(x_n))+c_2f(x_n+ah,y_n+bhf(x_n,y(x_n))))$$

这时，需要满足

$$\begin{cases} c_1+c_2=1 \\ 2c_2a=1 \\ 2c_2b=1 \\ \end{cases}$$

不难看出，这个方程组有无数多组解。

一个比较常用的二阶Runge-Kutta公式为

$$\begin{cases} y_{n+1}=y_n+\frac{h}{2}(k_1+k_2) \\ k_1=f(x_n,y_n) \\ k_2=f(x_n+h,y_n+hk_1) \end{cases}$$

一个比较常用的三阶Runge-Kutta公式为

$$\begin{cases} y_{n+1}=y_n+\frac{h}{6}(k_1+4k_2+k_3) \\ k_1=f(x_n,y_n) \\ k_2=f(x_n+\frac{h}{2},y_n+\frac{hk_1}{2}) \\ k_3=f(x_n+h,y_n-hk_1+2hk_2) \end{cases}$$

一个比较常用的四阶Runge-Kutta公式为

$$\begin{cases} y_{n+1}=y_n+\frac{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{hk_1}{2}) \\ k_3=f(x_n+\frac{h}{2},y_n+\frac{hk_2}{2}) \\ k_3=f(x_n+h,y_n+hk_3) \end{cases}$$

Runge-Kutta方法是一个在工程中应用很广泛的算法。

线性多步法

对于方程$\frac{\mathbf{d}y}{\mathbf{d}x}=f(x,y)$在区间$[x_{n-p},x_{n+1}]$上进行积分得

$$y(x_{n+1})=y(x_{n-p})+\int_{x_{n-p}}^{x_{n+1}}f(x,y)\mathrm{d}x$$

使用数值积分进行近似，从而构造线性多步格式。

如果使用积分节点$\lbrace x_n,x_{n-1}, \cdots ,x_{n-q} \rbrace$进行近似计算(详见插值型数值积分)，那么这就是一个显式求解格式。如果积分节点是$\lbrace x_{n+1},x_n,x_{n-1}, \cdots ,x_{n+1-q} \rbrace$，那么这就是一个隐式求解格式。

把$p=0$的格式称为Adams公式。

常微分方程组的数值解法

对于一个常微分方程组，可以把解常微分方程的各种方法平行地应用到方程组的求解中。

例如，对于下面一个方程组

$$\begin{cases} \frac{\mathbf{d}y}{\mathbf{d}x}=f(x,y,z) \\ \frac{\mathbf{d}z}{\mathbf{d}x}=g(x,y,z) \\ y(a)=y_0 \\ z(a)=z_0 \\ \end{cases}$$

其改进的欧拉公式为

$$\begin{pmatrix} \bar{y}_{n+1} \\ \bar{z}_{n+1} \\ \end{pmatrix} = \begin{pmatrix} y_n \\ z_n \end{pmatrix} +h \begin{pmatrix} f(x_n,y_n,z_n) \\ g(x_n,y_n,z_n) \\ \end{pmatrix} \\ \begin{pmatrix} y_{n+1} \\ z_{n+1} \end{pmatrix} = \begin{pmatrix} y_n \\ z_n \end{pmatrix} +\frac{h}{2}\left[ \begin{pmatrix} f(x_n,y_n,z_n) \\ g(x_n,y_n,z_n) \\ \end{pmatrix} + \begin{pmatrix} f(x_{n+1},\bar{y}_{n+1},\bar{z}_{n+1}) \\ g(x_{n+1},\bar{y}_{n+1},\bar{z}_{n+1}) \\ \end{pmatrix} \right]$$

高阶常微分方程的数值方法

对于高阶常微分方程，可以将其转换为一阶常微分方程组。

考虑一个二阶常微分方程

$$\begin{cases} \frac{\mathbf{d}^2y}{\mathbf{d}x^2}=f(x,y,y')\\ y(a)=\eta^{(0)} \\ y'(a)=\eta^{(1)} \\ \end{cases}$$

可以将其转换为

$$\begin{cases} \frac{\mathbf{d}y_1}{\mathbf{d}x}=y_2(x)\\ \frac{\mathbf{d}y_2}{\mathbf{d}x}=f(x,y_1(x),y_2(x)) \\ y_1(a)=\eta^{(0)} \\ y_2(a)=\eta^{(1)} \\ \end{cases}$$

对于更高阶的方程，也可以进行类似的转换。