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)\)

例题

多步法例子

其他

\[\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 法、梯形法均无条件稳定

例题

*本章做题的一般流程

此次考试中，本章内容较难且多，大部分难以理解，且作为最后一部分内容，从各方面来说对于首次接触的同学来说都是不太友好的。因此，老师贴心地对 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\)，即