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
多步法例子
其他
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\)