2022-2023 春学期 矩阵与数值分析 C7 常微分方程的数值解法

C7 常微分方程的数值解法

问题描述

一阶常微分方程的初值问题

$$\left\{ \begin{array}{l} u'=f(t,u),\;a\leq t\leq b\\ u(a)=u_0 \end{array} \right.$$

解的存在唯一性：若 $f(t,u)$ 满足 Lipschitz 条件，即存在常数 L，对任意 $t\in[a,b]$ 均有 $|f(t,u)-f(t,\bar{u})|\leq L|u-\bar{u}|$ 则问题的解存在且唯一。

但初值问题大多数情况下无法求得解析解，因此求其数值解；使用离散化的方法，在一系列预先取定的离散点（节点）中求出未知函数的近似值作为初值问题的数值解。

7.1 基于数值积分的数值方法

思想

将节点取为 $t_n=a+nh\quad h=\frac{b-a}{N},\quad n=0,1,2,\cdots,N$，每个节点区间求积分

$$\left\{ \begin{array}{l} u(t_n+1)=u(t_n)+\int^{t_{n+1}}_{t_n}f(t,u(t))dt\\ u(t_0)=u_0 \end{array} \right.$$

如果 $y(t_n)$ 的近似值 $u_n$ 已经求出，则计算右端项的数值积分可求出 $u(t_{n+1})$ 的近似值 $u_{n+1}$

Euler 法

对有段积分项使用左矩形公式，则得

$$u(t_{n+1})=u(t_n)+\int_{t_n}^{t_{n+1}}f(t,u(t))dt\approx u(t_n)+hf(t_n,u(t_n))$$

令

$$u_{n+1}=u_n+hf(t_n,u_n)\qquad n=0,1,2,\cdots,N-1$$

上式称为 Euler 求解公式，又称矩形公式。

该公式具有一阶精度，局部截断误差主项为 $1/2$

隐 Euler 法

对右端积分项使用右矩形求积公式，则得

$$u(t_{n+1})=u(t_n)+\int_{t_n}^{t_{n+1}}f(t,u(t))dt\approx u(t_n)+hf(t_{n+1},u(t_{n+1}))$$

令

$$u_{n+1}=u_n+hf(t_{n+1},u_{n+1})\qquad n=0,1,2,\cdots,N-1$$

上式称为隐 Euler 公式，又称右矩形公式，或向后 Euler 公式

隐式公式需要求解方程，或利用迭代法求解；可通过移项的方式将其显式化

该公式具有一阶精度，局部截断误差主项为 $-1/2$

梯形法

对右端积分项使用梯形求积公式，则得

$$u(t_{n+1})=u(t_n)+\int_{t_n}^{t_{n+1}}f(t,u(t))dt\approx u(t_n)+\frac{h}{2}[f(t_{n},u(t_{n}))+f(t_{n+1},u(t_{n+1}))]$$

令

$$u_{n+1}=u_n+\frac{h}{2}[f(t_{n},u(t_{n}))+f(t_{n+1},u(t_{n+1}))]\qquad n=0,1,2,\cdots,N$$

上式称为梯形公式，简称梯形法，梯形公式可看作 Euler 公式与隐式 Euler 公式的算数平均

该公式具有二阶精度，局部截断误差主项为 $-1/12$

改进的 Euler 法

为避免求解非线性代数方程，可以用 Euler 法将其显化，建立预测——校正系统：

$$\left\{ \begin{array}{l} u(t_0)=u_0\\ \bar{u}_{n+1}=u_n+hf(t_n,u_n)\\ u_{n+1}=u_n+\frac{h}{2}(f(t_n,u_n)+f(t_{n+1},\bar{u}_{n+1})) \end{array} \right.$$

求解公式称为改进的 Euler 法，其中 称为预测值， 称为校正值

其求解顺序为：

$$u_0\to \bar{u}_1\to u_1\to \bar{u}_2 \to u_2\to\cdots\to\bar{u}_N\to u_N$$

改进的 Euler 法还可写为：

$$u_{n+1}=u_n+\frac{h}{2}[f(t_n,u_n)+f(t_{n+1},u_n+hf(t_n,u_n) )]$$

该公式具有二阶精度

截断误差与精度

局部截断误差：假设 $u_i=u(t_i),\;i=0,1,2,\cdots,n-1$ 则称 $R_n(h)=u(t_n)-u_n$ 为求解公式第 n 步的局部截断误差

整体截断误差：$E_n(h)=u(t_n)-u_n=\sum\limits^n_{t=1}R_t(h)$ 为求解公式在 $t_n$ 点上的整体截断误差

求解公式具有 p 阶精度：

• 求解公式的局部截断误差：$R_n(h)=O(h^{p+1})$
• 求解公式的整体截断误差：$E_n(h)=O(h^p)$

例题

1

多步法例子

其他

$$\frac{1}{2}\times(-50)h<1$$

7.2 基于 Taylor 展式的数值方法

Runge-Kutta 法（不考

待定系数法

线性多步法（k 步线性法）的一般公式：

$$\sum\limits^k_{j=0}\alpha_ju_{n+1}=h\sum\limits^k_{j=0}\beta_jf_{n+j},\;\alpha_k\neq0$$

其中 $f_{n+j}=f(t_{n+j},u_{n+j})$, $\alpha_j,\beta_j$ 是常数，$\alpha_0$ 和 $\beta_0$ 不同时为 0

定理（多步法性质定理）：多步法和下列三个性质等价：

• p阶精度：

  $$c_0=c_1=c_2\cdots=c_p=0$$

• *对每个次数 $\leq p$ 的多项式，$L[f_p(t);h]=0$

• *对一切 $u(t)\in C^{p+1},L[u(t);h]=O(h^{p+1})$

局部截断误差主项：$c_{p+1}h^{p+1}u^{(p+1)}(t)$ 称为局部截断误差主项；其中，$c_{p+1}$ 称为局部截断误差主项系数，其整体截断误差 $E_n=O(h^p)$，故称为 p 阶 k 步法。

$$\left\{ \begin{array}{l} c_0=\alpha_0+\alpha_1+\cdots+\alpha_k=0\\ c_1=\alpha_1+2\alpha_2+\cdots+k\alpha_k-(\beta_0+\beta_1+\cdots+\beta_k)=0\\ \qquad\qquad\qquad\vdots\\ c_p=\frac{1}{p!}(\alpha_1+2^p\alpha_2+\cdots+k^p\alpha_k)-\frac{1}{(p-1)!}(\beta_1+2\beta^{p-1}\beta_2+\cdots+k^{p-1}\beta_k)=0\\ c_{p+1}\neq0 \end{array} \right.$$

有 p+1 个方程，2k+1 个未知量，因此 $p\leq 2k$；线性 k 步法最高可达到 2k 阶精度（整体截断误差）

观察可知，$c_i$ 中，$\beta$ 部分的系数均与 $c_{i-1}$ 中 $\alpha$部分的系数一致。

收敛性

对于单步法，当方法的阶 $p\geq1$ 时，有整体误差

$$E_n=u(t_n)-u_n=O(h^p)$$

故有 $\lim\limits_{h\to0}E_n=0$，因此方法是收敛的

对多步法，若方法是 k 步 p 阶法，

$$\sum\limits_{j=0}^k\alpha_ju_{n+j}=h\sum\limits^k_{j=0}\beta_jf_{n+j},\;\alpha_k\neq0$$

则多步法的第一特征多项式为：

$$\rho(\lambda)=\sum\limits_{j=0}^k\alpha_j\lambda^j$$

第二特征多项式为：

$$\sigma(\lambda)=\sum\limits_{j=0}^k\beta_j\lambda^j$$

若多步法的第一特征多项式 $\rho(\lambda)$ 的所有根（复根）在单位圆或圈上（$|\lambda|\leq1$），且位于单位圆周上的根都是单根，则称多步法满足根条件。若线性多步法的阶 $p\geq1$，且满足根条件，则方法是收敛的。

所以，收敛需要同时满足如下两个条件：

• 依据第一特征多项式 $\rho(\lambda)=\sum\limits_{j=0}^k\alpha_j\lambda^j=0$，解出 $\lambda$，判断是否有 $|\lambda|\leq1$
• 判断该多步法的阶是否有：$p\geq1$

绝对稳定区间

定义：一个数值方法用于求解模型问题 $u’=\mu u$，若在平面中的某一区域 D 中方法都是绝对稳定的，而在区域 D 外，方法是不稳定的，则称 D 是方法的绝对稳定区域，它与实轴的交称为绝对稳定区间

特征方程：

$$\rho(\lambda)-\bar{h}\sigma(\lambda)=0$$

其中

$$\rho(\lambda)=\sum\limits^k_{j=0}\alpha_j\lambda^j\qquad \sigma(\lambda)=\sum\limits^k_{j=0}\beta_j\lambda^j\qquad \bar{h}=\mu h$$

由于 $\mu<0$，所以 $\bar{h}<0$

绝对稳定域：若特征方程 $\rho(\lambda)-\bar{h}\sigma(\lambda)=0$ 的根都在单位圆内（$|\lambda|<1$），则线性多步法关于 $\bar{h}=\mu h$ 绝对稳定，其绝对稳定域是复平面 $\bar{h}$ 上的区域：

$$D=\{\bar{h}|\lambda_j(\bar{h})<1,\quad j=1,2,\cdots,k\}$$

绝对稳定域的使用判别法：实系数二次方程 $\lambda^2-b\lambda-c=0$ 的根在单位圆的充分必要条件为：

$$|b|<1-c<2$$

可通过上述方法得到：隐 Euler 法、梯形法 均无条件稳定

例题

1

2

3

4

5

6

7

*本章做题的一般流程

此次考试中，本章内容较难且多，大部分难以理解，且作为最后一部分内容，从各方面来说对于首次接触的同学来说都是不太友好的。因此，老师贴心地对 Runge-Kutta 部分不做考核要求。

此处总结部分大题中可能涉及到的套路，小题概念部分可参考课本或课程 PPT。

首先，将整理所给条件，分别将 $u,f$ 移动到等号两边，成为下式的形式

$$\sum\limits_{j=0}^k\alpha_ju_{n+j}=h\sum\limits^k_{j=0}\beta_jf_{n+j}$$

• 求方法阶数、截断误差主项系数、主项

  构造线性方程组，得到有 $c_0=c_1=\cdots=c_p=0$ ，而 $c_{p+1}\neq0$ 为局部截断误差主项的系数，具体如下：

  $$\left\{ \begin{array}{l} c_0=\alpha_0+\alpha_1+\cdots+\alpha_k=0\\ c_1=\alpha_1+2\alpha_2+\cdots+k\alpha_k-(\beta_0+\beta_1+\cdots+\beta_k)=0\\ \qquad\qquad\qquad\vdots\\ c_p=\frac{1}{p!}(\alpha_1+2^p\alpha_2+\cdots+k^p\alpha_k)-\frac{1}{(p-1)!}(\beta_1+2\beta^{p-1}\beta_2+\cdots+k^{p-1}\beta_k)=0\\ c_{p+1}\neq0 \end{array} \right.$$

  得到局部误差主项为 $c^{p+1}h^{p+1}u^{p+1}(t_n)$

  此方法为 k 步 p 阶法

• 讨论收敛性

  依据第一特征多项式 $\rho(\lambda)=\sum\limits_{j=0}^k\alpha_j\lambda^j=0$，解出 $\lambda$，判断是否有 $|\lambda|\leq1$；并判断是否有阶数 $p\geq1$；若上述条件均满足，则收敛

• 求绝对稳定区间

  构造特征方程

  $$\rho(\lambda)-\bar{h}\sigma(\lambda)=0$$

  并转换为实系数二次方程：

  $$\lambda^2-b\lambda-c=0$$

  求解不等式

  $$|b|<1-c<2$$

  得到 $\bar{h}$ 的范围即为绝对稳定区间，需要注意的是 $\bar{h}<0$ 被认为是必须满足的已知条件

可结合上节例7或本节下述例题部分理解该过程

例题

证明求解一阶常微分方程初值问题：$u’=f(t,u),u(0)=u_0$ 的差分格式 $u_{n+2}-u_{n+1}=\frac{h}{12}(5f_{n+2}+8f_{n+1}-f_n)$ 收敛并求其局部截断误差主项，绝对稳定区间

解：

依题意可得

$$\alpha_0=0\quad\alpha_1=-1\quad\alpha_2=1$$ $$\beta_0=\frac{-1}{12}\quad\beta_1=\frac{8}{12}\quad\beta_2=\frac{5}{12}$$

从而可以构造线性方程组，得到

$$\left\{ \begin{array}{l} c_0=1-1=0\\ c_1=(-1+2\times1)-(\frac{-1}{12}+\frac{8}{12}+\frac{5}{12})=0\\ c_2=\frac{1}{2}(-1+1*2^2)-(\frac{8}{12}+2\times\frac{5}{12})=0\\ c_3=\frac{1}{6}(-1+1*2^3)-\frac{1}{2}(\frac{8}{12}+2^2\times\frac{5}{12})=0\\ c_4=\frac{1}{24}(-1+1*2^4)-\frac{1}{6}(\frac{8}{12}+2^3\times\frac{5}{12})=\frac{-1}{24}\neq0 \end{array} \right.$$

因此，此方法为隐式二步三阶法，其局部截断误差主项为 $\frac{-1}{24}h^4u^{(4)}(t_n)$

由差分格式可构造第一特征多项式

$$\rho(\lambda)=\lambda^2-\lambda$$

第二特征多项式

$$\sigma(\lambda)=\frac{5}{12}\lambda^2+\frac{8}{12}\lambda-\frac{1}{12}$$

令 $\rho(\lambda)=\lambda^2-\lambda=0$

得 $\lambda_1=0,\;\lambda_2=1$，则其特征值满足根条件 $|\lambda|<1$，且阶数 $p\geq1$，因此该方法收敛

构造特征方程：

$$\rho(\lambda)-\bar{h}\sigma(\lambda)= (1-\frac{5}{12}\bar{h})\lambda^2-(1+\frac{8}{12}\bar{h})\lambda+\frac{1}{12}\bar{h} =0$$

将特征方程二次项系数化为 1，可得

$$\lambda^2- (\frac{12+8\bar{h}}{12-5\bar{h}})\lambda- \frac{(-\bar{h})}{12-5\bar{h}}=0$$

求解不等式 $|b|<1-c<2$，即

$$\left| (\frac{12+8\bar{h}}{12-5\bar{h}}) \right|< 1-\frac{(-\bar{h})}{12-5\bar{h}}= \frac{12-4\bar{h}}{12-5\bar{h}}<2$$

注意到已知 $\bar{h}<0$

而 $1-\frac{(-\bar{h})}{12-5\bar{h}}=
\frac{12-4\bar{h}}{12-5\bar{h}}<2$ 自然成立，

再求解 $
\left|
(\frac{12+8\bar{h}}{12-5\bar{h}})
\right|<\frac{12-4\bar{h}}{12-5\bar{h}}$

得到 $-12+4\bar{h}<12+8\bar{h}<12-4\bar{h}$，

而 $12+8\bar{h}<12-4\bar{h}$ 自然成立，

即有 $-3+\bar{h}<3+2\bar{h}$，

可得其绝对稳定区间为 $-6<\bar{h}<0$

END