ch09+ch01+常微分方程初、边值问题数值解法+硕博讲义
第1章 常微分方程初、边值问题数值解法 1.1 引言 详细讲解
作为拥有多年教学经验的资深高级研究员，我将从历史脉络、理论基础、核心问题、方法框架四个维度，对本节内容进行系统、深入的讲解，并给出关键定理的完整证明。
一、微分方程的历史起源与早期里程碑式研究
微分方程与微积分是孪生学科，几乎同时诞生于17世纪下半叶。它的出现源于人类对自然界运动规律的定量描述需求，早期数学家的开创性工作奠定了整个学科的基础。
1. 艾萨克·牛顿（Isaac Newton, 1643-1727）的贡献
牛顿在1671年完成的《流数术与无穷级数》（1736年正式出版）中，首次系统地将微分方程作为研究运动的工具。他研究的第一个一阶微分方程是：
核心意义：
- 这是人类历史上第一个被系统研究的非线性一阶微分方程
- 牛顿首创了用级数解法和积分法求解微分方程的思想，为后续解析解法奠定了基础
- 该方程源于牛顿对天体运动和力学系统的研究，体现了微分方程与物理学的天然联系
2. 戈特弗里德·威廉·莱布尼茨（Gottfried Wilhelm Leibniz, 1646-1716）的贡献
莱布尼茨约于1676年在研究几何中的反切线问题时，建立了如下微分方程模型：
反切线问题的几何意义：
已知曲线在任意一点的切线斜率满足某种规律，求该曲线的方程。这是微分方程几何应用的经典问题。
求解过程：
这是一个可分离变量的微分方程，分离变量得：
两边积分：
令\(y = a\sin\theta\)，则\(dy = a\cos\theta d\theta\)，代入得：
代回\(\theta = \arcsin\frac{y}{a}\)，得到通解：
核心意义：
- 莱布尼茨创立了现代微分方程的符号体系（如\(dy/dx\)），极大地推动了学科的发展
- 他系统地研究了可分离变量方程和齐次方程的解法，建立了初等积分法的基本框架
3. 莱昂哈德·欧拉（Leonhard Euler, 1707-1783）的贡献
欧拉于1744年在《寻求具有某种极大或极小性质的曲线的方法》中，借助二阶微分方程：
给出了变分法中极小问题：
的一般解（欧拉-拉格朗日方程）。
核心意义：
- 欧拉将微分方程的应用从力学和几何扩展到了变分法领域
- 他创立了欧拉法，这是第一个系统的常微分方程数值解法，为数值分析学科的诞生奠定了基础
- 欧拉还系统地研究了线性微分方程的理论，建立了线性微分方程解的结构定理
4. 亚历克西斯·克劳德·克莱罗（Alexis Claude Clairaut, 1713-1765）的贡献
克莱罗于1734年在研究长方形框的移动问题时，建立了第一个隐式微分方程：
解的结构：
该方程有两类解：
- 通解（直线族）：\(y = Cx - f(C)\)，其中\(C\)为任意常数
- 奇解（包络曲线）：通过消去参数\(C\)，从方程组\(\begin{cases}y = Cx - f(C) \\ 0 = x - f'(C)\end{cases}\)得到
例子：当\(f(y') = (y')^2\)时，克莱罗方程为\(y - xy' + (y')^2 = 0\)
- 通解：\(y = Cx - C^2\)（直线族）
- 奇解：由\(0 = x - 2C\)得\(C = x/2\)，代入通解得\(y = x^2/4\)（抛物线，是直线族的包络）
核心意义： - 克莱罗首次发现了微分方程的奇解现象，揭示了微分方程解的复杂性
- 隐式微分方程的研究推动了微分方程几何理论的发展
二、解析解与渐近分析的局限性
引言中指出："在这些方程中，仅有很少的一部分能通过初等积分法给出其通解或通积分。"这是微分方程理论中一个极其重要的结论，也是数值解法产生的根本原因。
1. 初等积分法的局限性
初等函数：由常数、基本初等函数（幂函数、指数函数、对数函数、三角函数、反三角函数）经过有限次四则运算和复合运算得到的函数。
初等积分法：通过积分运算将微分方程的解表示为初等函数的方法。
为什么绝大多数微分方程没有初等解？
- 1841年，法国数学家约瑟夫·刘维尔（Joseph Liouville）证明了：黎卡提方程\(y' = P(x)y^2 + Q(x)y + R(x)\)，除非能找到一个特解，否则一般没有初等解。
- 即使是形式非常简单的微分方程，如\(y' = x^2 + y^2\)，也没有初等解。
- 随着微分方程阶数的提高和非线性程度的增强，存在初等解的方程比例急剧下降。
例子：著名的范德波尔方程\(y'' + \mu(y^2 - 1)y' + y = 0\)（\(\mu > 0\)），描述了非线性电子振荡电路的行为，至今没有找到初等解。
2. 渐近分析的局限性
渐近分析是研究微分方程解在极限情况下（如\(x \to \infty\)或\(x \to 0\)）行为的方法。它在理论研究中有重要作用，但存在明显的局限性：
- 只能给出解在局部区域的近似行为，不能给出全域的精确解
- 对于非极限区域的解，渐近分析的误差可能非常大
- 许多实际问题需要知道解在整个区间上的具体数值，而不仅仅是极限行为
三、数值解法的兴起与重要性
在计算机迅猛发展的今天，微分方程的数值求解已经成为科学与工程计算中最重要的工具之一。
1. 数值解法的核心优势
- 通用性：可以求解绝大多数无法用解析方法和渐近方法解决的微分方程
- 精确性：通过减小步长或使用高阶方法，可以将误差控制在任意需要的范围内
- 高效性：借助现代计算机，可以在短时间内求解大规模、复杂的微分方程组
2. 数值解法与理论分析的关系
- 数值模拟需要理论检验：数值解法得到的结果必须通过理论分析或实验验证，以确保其正确性
- 理论分析可以通过数值方法简化：许多复杂的理论问题可以通过数值模拟得到直观的理解，从而推动理论的发展
四、本章核心问题与理论基础
本章主要研究一阶常微分方程初值问题的数值解法，其数学模型为：
其中\(f(x, y)\)是定义在区域\(D = [a, b] \times \mathbb{R}\)上的二元函数，\(y_a\)是给定的初始值。
1. 解的存在唯一性定理（皮卡-林德勒夫定理）
这是所有数值解法的理论基础。如果初值问题的解不存在或者不唯一，那么数值解法就失去了意义。
定理条件：
如果函数\(f(x, y)\)满足：
- 连续性：\(f(x, y)\)在闭矩形区域\(R = \{(x, y) \mid |x - a| \leq A, |y - y_a| \leq B\}\)上连续
- 利普希茨条件：\(f(x, y)\)关于\(y\)满足利普希茨条件，即存在常数\(L > 0\)，使得对任意\((x, y_1), (x, y_2) \in R\)，有\[|f(x, y_1) - f(x, y_2)| \leq L|y_1 - y_2| \]
其中\(L\)称为利普希茨常数。
定理结论：
初值问题(1.1.1)-(1.1.2)在区间\(|x - a| \leq h\)上存在唯一的连续解\(y = y(x)\)，其中
完整证明过程：
我们将证明分为五个步骤：
步骤1：将初值问题转化为等价的积分方程
如果\(y = y(x)\)是初值问题(1.1.1)-(1.1.2)的解，那么对(1.1.1)两边从\(a\)到\(x\)积分，并利用初始条件(1.1.2)，得：
反之，如果\(y = y(x)\)是积分方程(1.1.3)的连续解，那么对(1.1.3)两边求导，得：
且\(y(a) = y_a + \int_a^a f(t, y(t)) dt = y_a\)，即\(y = y(x)\)是初值问题(1.1.1)-(1.1.2)的解。
因此，初值问题(1.1.1)-(1.1.2)与积分方程(1.1.3)是等价的。
步骤2：构造皮卡迭代序列
我们构造如下的迭代序列：
这个序列称为皮卡迭代序列。
步骤3：证明皮卡迭代序列在区间\(|x - a| \leq h\)上一致收敛
首先，我们用数学归纳法证明：对所有\(n \geq 0\)，\(y_n(x)\)在区间\(|x - a| \leq h\)上连续，且满足\(|y_n(x) - y_a| \leq B\)。
- 当\(n = 0\)时，\(y_0(x) = y_a\)显然连续，且\(|y_0(x) - y_a| = 0 \leq B\)。
- 假设当\(n = k\)时，结论成立，即\(y_k(x)\)在区间\(|x - a| \leq h\)上连续，且\(|y_k(x) - y_a| \leq B\)。
那么当\(n = k+1\)时，\(y_{k+1}(x) = y_a + \int_a^x f(t, y_k(t)) dt\)。
由于\(f(t, y_k(t))\)是连续函数的复合，因此是连续的，其积分也是连续的，故\(y_{k+1}(x)\)连续。
又因为\[|y_{k+1}(x) - y_a| = \left|\int_a^x f(t, y_k(t)) dt\right| \leq \int_a^x |f(t, y_k(t))| dt \leq M|x - a| \leq Mh \leq B \]所以结论对\(n = k+1\)也成立。
由数学归纳法，结论对所有\(n \geq 0\)成立。
接下来，我们考虑级数：
其部分和为\(S_N(x) = y_0(x) + \sum_{n=0}^{N-1} [y_{n+1}(x) - y_n(x)] = y_N(x)\)。
因此，证明皮卡迭代序列一致收敛等价于证明级数(1.1.5)一致收敛。
我们用数学归纳法证明：对所有\(n \geq 0\)，有
- 当\(n = 0\)时，\[|y_1(x) - y_0(x)| = \left|\int_a^x f(t, y_0(t)) dt\right| \leq \int_a^x |f(t, y_a)| dt \leq M|x - a| \]结论成立。
- 假设当\(n = k\)时，结论成立，即\[|y_k(x) - y_{k-1}(x)| \leq \frac{M L^{k-1}}{k!} |x - a|^k \]那么当\(n = k+1\)时，\[\begin{align*} |y_{k+1}(x) - y_k(x)| &= \left|\int_a^x [f(t, y_k(t)) - f(t, y_{k-1}(t))] dt\right| \\ &\leq \int_a^x |f(t, y_k(t)) - f(t, y_{k-1}(t))| dt \\ &\leq L \int_a^x |y_k(t) - y_{k-1}(t)| dt \\ &\leq L \int_a^x \frac{M L^{k-1}}{k!} |t - a|^k dt \\ &= \frac{M L^k}{k!} \cdot \frac{|x - a|^{k+1}}{k+1} \\ &= \frac{M L^k}{(k+1)!} |x - a|^{k+1} \end{align*}\]结论成立。
由数学归纳法，不等式(1.1.6)对所有\(n \geq 0\)成立。
因为\(|x - a| \leq h\)，所以
而级数\(\sum_{n=0}^{\infty} \frac{M L^n h^{n+1}}{(n+1)!} = \frac{M}{L} (e^{Lh} - 1)\)是收敛的正项级数。
根据魏尔斯特拉斯判别法，级数(1.1.5)在区间\(|x - a| \leq h\)上一致收敛。
因此，皮卡迭代序列\(\{y_n(x)\}\)在区间\(|x - a| \leq h\)上一致收敛，设其极限为\(y(x)\)，即
步骤4：证明极限函数\(y(x)\)是积分方程(1.1.3)的解
由于\(\{y_n(x)\}\)一致收敛于\(y(x)\)，且\(f(x, y)\)关于\(y\)满足利普希茨条件，因此\(f(x, y_n(x))\)一致收敛于\(f(x, y(x))\)。
对迭代公式(1.1.4)两边取极限，得：
即\(y(x)\)是积分方程(1.1.3)的解，从而也是初值问题(1.1.1)-(1.1.2)的解。
步骤5：证明解的唯一性
假设\(y(x)\)和\(z(x)\)都是初值问题(1.1.1)-(1.1.2)在区间\(|x - a| \leq h\)上的解，那么它们都满足积分方程：
两式相减，得：
两边取绝对值，得：
令\(u(x) = \int_a^x |y(t) - z(t)| dt\)，则\(u'(x) = |y(x) - z(x)|\)，且\(u(a) = 0\)。
不等式(1.1.7)变为：
即
两边乘以\(e^{-Lx}\)，得：
因此，函数\(u(x) e^{-Lx}\)在区间\(|x - a| \leq h\)上是单调递减的。
所以，对任意\(x \geq a\)，有
即\(u(x) \leq 0\)。
但\(u(x) = \int_a^x |y(t) - z(t)| dt \geq 0\)，因此\(u(x) = 0\)对所有\(x \in [a, a+h]\)成立。
从而\(|y(x) - z(x)| = u'(x) = 0\)对所有\(x \in [a, a+h]\)成立，即\(y(x) = z(x)\)。
同理可证，对\(x \in [a-h, a]\)，也有\(y(x) = z(x)\)。
因此，初值问题(1.1.1)-(1.1.2)的解是唯一的。
定理的重要说明：
- 利普希茨条件是保证解唯一的充分条件，但不是必要条件。例如，方程\(y' = \sqrt{y}\)，\(y(0) = 0\)，有两个解\(y = 0\)和\(y = x^2/4\)，但\(f(x, y) = \sqrt{y}\)在\(y = 0\)处不满足利普希茨条件。
- 如果\(f(x, y)\)在区域\(D\)上关于\(y\)的偏导数\(f_y(x, y)\)存在且有界，那么\(f(x, y)\)关于\(y\)满足利普希茨条件，且利普希茨常数\(L = \max_{(x, y) \in D} |f_y(x, y)|\)。这是验证利普希茨条件最常用的方法。
2. 数值解法的基本思想：离散化
数值解法的核心思想是将连续的问题转化为离散的问题。具体来说：
- 区间离散化：将求解区间\([a, b]\)分成\(n\)个相等的子区间，分点为：\[x_0 = a, \quad x_1 = a + h, \quad x_2 = a + 2h, \quad \dots, \quad x_n = a + nh = b \]其中\(h = (b - a)/n\)称为步长。
- 函数离散化：我们的目标不是求连续函数\(y = y(x)\)，而是求它在离散点\(x_0, x_1, \dots, x_n\)上的近似值\(y_0, y_1, \dots, y_n\)，其中\(y_0 = y_a\)是已知的初始值。
- 递推公式：通过某种数值方法，建立从\(y_i\)计算\(y_{i+1}\)的递推公式，然后从\(y_0\)出发，依次计算\(y_1, y_2, \dots, y_n\)。
数值解法的分类：
- 单步法：计算\(y_{i+1}\)时，只需要用到前一个点的信息\(y_i\)，如欧拉法、龙格-库塔法。
- 多步法：计算\(y_{i+1}\)时，需要用到前多个点的信息\(y_i, y_{i-1}, \dots, y_{i-k+1}\)，如亚当斯法、米尔恩法。
- 显式方法：递推公式中\(y_{i+1}\)可以直接由已知量计算得到。
- 隐式方法：递推公式中\(y_{i+1}\)出现在等式两边，需要通过解方程才能得到。
五、本章内容框架
本章将按照以下顺序系统地介绍常微分方程的数值解法：
- 一阶常微分方程初值问题的基本数值方法：
- 欧拉法（显式欧拉法、隐式欧拉法、梯形法）
- 龙格-库塔法（二阶龙格-库塔法、四阶龙格-库塔法）
- 线性多步法（亚当斯显式法、亚当斯隐式法、预测-校正法）
- 数值方法的理论分析：
- 截断误差与阶
- 收敛性与稳定性
- 绝对稳定性与绝对稳定区域
- 高阶常微分方程与一阶微分方程组的数值解法：
- 高阶方程转化为一阶方程组
- 一阶方程组的数值解法
- 刚性常微分方程（组）的数值解法：
- 刚性问题的定义与特点
- 刚性稳定方法（如向后差分法BDF）
- 常微分方程边值问题的数值解法：
- 打靶法
- 有限差分法
- Hamilton系统的辛几何算法：
- Hamilton系统的基本性质
- 辛几何算法的基本思想与构造
六、核心知识点归纳总结
为了便于大家系统掌握本节内容，我将核心知识点归纳为以下表格：
| 知识点类别 | 具体内容 | 核心要点 |
|---|---|---|
| 历史里程碑 | 牛顿(1671) | 研究第一个一阶微分方程，首创级数解法和积分法 |
| 莱布尼茨(1676) | 创立现代微分方程符号体系，研究反切线问题，建立可分离变量方程解法 | |
| 欧拉(1744) | 建立欧拉-拉格朗日方程，创立欧拉法，奠定数值分析基础 | |
| 克莱罗(1734) | 发现隐式微分方程的奇解现象，推动微分方程几何理论发展 | |
| 核心问题 | 一阶常微分方程初值问题 | 数学模型：\(\frac{dy}{dx}=f(x,y), y(a)=y_a, x\in[a,b]\) |
| 理论基础 | 解的存在唯一性定理(皮卡-林德勒夫定理) | 条件：\(f(x,y)\)连续且关于\(y\)满足利普希茨条件 结论：在区间\(|x-a|\leq h\)上存在唯一解 证明方法：皮卡迭代法 |
| 解析解局限性 | 初等积分法 | 绝大多数微分方程没有初等解，如黎卡提方程\(y'=P(x)y^2+Q(x)y+R(x)\) |
| 渐近分析 | 只能给出局部极限行为，不能给出全域精确解 | |
| 数值解法基本思想 | 离散化 | 1. 区间离散化：将\([a,b]\)分为\(n\)个等距子区间 2. 函数离散化：求\(y(x)\)在离散点\(x_i\)上的近似值\(y_i\) 3. 递推公式：建立从\(y_i\)计算\(y_{i+1}\)的公式 |
| 数值解法分类 | 按步数分 | 单步法(欧拉法、龙格-库塔法)、多步法(亚当斯法) |
| 按显隐式分 | 显式方法、隐式方法 | |
| 本章内容框架 | 初值问题基本方法 | 欧拉法、龙格-库塔法、线性多步法 |
| 理论分析 | 截断误差、收敛性、稳定性 | |
| 扩展问题 | 高阶方程与方程组、刚性问题、边值问题、辛几何算法 |
七、学习建议
- 深刻理解解的存在唯一性定理：这是所有数值解法的理论基础，必须掌握其条件、结论和证明过程。
- 掌握离散化的基本思想：数值解法的本质是将连续问题转化为离散问题，理解这一点是学习各种数值方法的关键。
- 注重理论与实践结合：在学习数值方法的同时，要通过编程实现这些方法，加深对算法原理和误差分析的理解。
- 关注方法的适用范围：不同的数值方法有不同的优缺点和适用范围，要根据具体问题选择合适的方法。
1.1 节后续内容 详细讲解与证明
作为拥有60多年教学经验的资深高级研究员，我将对教材中这部分内容进行系统、深入的讲解，补充教材中省略的所有关键证明，并对核心概念进行深度解读。
一、存在唯一性定理（教材版本）的深度解读
教材给出的定理1.1.1是皮卡-林德勒夫定理的全局版本，与我上一讲中介绍的局部版本相比，条件更强，但结论也更具一般性。
定理原文回顾
定理1.1.1（存在唯一解） 如果方程(1.1.1)中的右端函数\(f(x,y)\)满足：
(i) \(f(x,y)\)是实值函数，
(ii) 函数\(f(x,y)\)在矩形区域\(\Omega = \{(x,y) \mid x \in [a,b], \ y \in \mathbb{R}\}\)内连续，
(iii) \(f(x,y)\)关于\(y\)满足Lipschitz条件：存在正常数\(L\)，使得对任意\(x \in [a,b]\)，均成立不等式
则问题(1.1.1)-(1.1.2)存在唯一的解\(y(x) \in C^1[a,b]\)。
与局部版本的对比
| 版本 | 区域 | 结论 | 适用范围 |
|---|---|---|---|
| 局部版本（上一讲） | 有限闭矩形\(R = \{(x,y) \mid |x-a| \leq A, |y-y_a| \leq B\}\) | 在区间\(|x-a| \leq h\)上存在唯一解 | 所有满足局部Lipschitz条件的函数 |
| 全局版本（教材） | 无限带形区域\(\Omega = [a,b] \times \mathbb{R}\) | 在整个区间\([a,b]\)上存在唯一解 | 满足全局Lipschitz条件的函数 |
全局版本的完整证明
教材说"上述定理的证明可参阅有关微分方程理论的教材"，这里我给出完整、严谨的证明。
证明：我们将证明分为两个部分：存在性和唯一性。
1. 唯一性证明
假设\(y(x)\)和\(z(x)\)都是初值问题在区间\([a,b]\)上的解，那么它们都满足积分方程：
两式相减得：
两边取绝对值，并利用Lipschitz条件：
令\(u(x) = |y(x) - z(x)|\)，则上式变为：
这是一个Gronwall不等式。我们来证明\(u(x) \equiv 0\)。
令\(v(x) = \int_a^x u(t) dt\)，则\(v'(x) = u(x)\)，且\(v(a) = 0\)。不等式(1.1.8)变为：
即
两边乘以\(e^{-Lx}\)（积分因子），得：
因此，函数\(v(x) e^{-Lx}\)在区间\([a,b]\)上单调递减。所以对任意\(x \in [a,b]\)，有：
即\(v(x) \leq 0\)。但\(v(x) = \int_a^x u(t) dt \geq 0\)，故\(v(x) \equiv 0\)。
从而\(u(x) = v'(x) \equiv 0\)，即\(y(x) \equiv z(x)\)。唯一性得证。
2. 存在性证明
我们仍然使用皮卡迭代法，但需要证明迭代序列在整个区间\([a,b]\)上一致收敛。
构造皮卡迭代序列：
第一步：证明所有\(y_n(x)\)在\([a,b]\)上连续
- \(y_0(x) = y_a\)显然连续。
- 假设\(y_n(x)\)连续，则\(f(t, y_n(t))\)连续（因为\(f\)连续），其积分\(y_{n+1}(x)\)也连续。
由数学归纳法，所有\(y_n(x)\)在\([a,b]\)上连续。
第二步：证明迭代序列一致收敛
考虑级数：
其部分和为\(S_N(x) = y_N(x)\)。我们用数学归纳法证明：
其中\(M = \max_{x \in [a,b]} |f(x, y_a)|\)。
- 当\(n=0\)时：\[|y_1(x) - y_0(x)| = \left| \int_a^x f(t, y_a) dt \right| \leq M (x - a) \]结论成立。
- 假设当\(n=k\)时结论成立，即：\[|y_k(x) - y_{k-1}(x)| \leq \frac{M L^{k-1}}{k!} (x - a)^k \]那么当\(n=k+1\)时：\[\begin{align*} |y_{k+1}(x) - y_k(x)| &= \left| \int_a^x [f(t, y_k(t)) - f(t, y_{k-1}(t))] dt \right| \\ &\leq \int_a^x |f(t, y_k(t)) - f(t, y_{k-1}(t))| dt \\ &\leq L \int_a^x |y_k(t) - y_{k-1}(t)| dt \\ &\leq L \int_a^x \frac{M L^{k-1}}{k!} (t - a)^k dt \\ &= \frac{M L^k}{k!} \cdot \frac{(x - a)^{k+1}}{k+1} \\ &= \frac{M L^k}{(k+1)!} (x - a)^{k+1} \end{align*}\]结论成立。
由数学归纳法，不等式(1.1.9)对所有\(n \geq 0\)成立。
因为\(x \in [a,b]\)，所以\((x - a) \leq (b - a)\)，因此：
而级数\(\sum_{n=0}^{\infty} \frac{M L^n (b - a)^{n+1}}{(n+1)!} = \frac{M}{L} \left( e^{L(b-a)} - 1 \right)\)是收敛的正项级数。
根据魏尔斯特拉斯M判别法，级数\(\sum_{n=0}^{\infty} [y_{n+1}(x) - y_n(x)]\)在区间\([a,b]\)上一致收敛。因此，皮卡迭代序列\(\{y_n(x)\}\)在\([a,b]\)上一致收敛，设其极限为\(y(x)\)。
第三步：证明极限函数\(y(x)\)是解
由于\(\{y_n(x)\}\)一致收敛于\(y(x)\)，且\(f(x,y)\)关于\(y\)满足Lipschitz条件，因此\(f(x, y_n(x))\)一致收敛于\(f(x, y(x))\)。
对迭代公式两边取极限：
即\(y(x)\)是积分方程的解，从而也是初值问题的解。
第四步：证明\(y(x) \in C^1[a,b]\)
因为\(y(x)\)是积分方程的解，所以\(y(x)\)是连续函数的积分，故\(y(x)\)连续。又因为\(f(x, y(x))\)连续，所以\(y'(x) = f(x, y(x))\)连续，即\(y(x) \in C^1[a,b]\)。
存在性得证。
重要说明
- 全局Lipschitz条件是一个很强的条件，许多实际问题中的函数不满足这个条件。例如，\(f(x,y) = y^2\)在\(\mathbb{R}\)上不满足全局Lipschitz条件，因为其导数\(f_y = 2y\)是无界的。
- 对于不满足全局Lipschitz条件的函数，我们只能保证解在某个局部区间上存在，解可能在有限时间内"爆破"（趋向于无穷大）。例如，初值问题\(y' = y^2, y(0) = 1\)的解是\(y(x) = 1/(1-x)\)，它在\(x=1\)处爆破。
二、适定性（连续依赖性）的深度讲解
适定性是微分方程理论中一个极其重要的概念，它是数值解法可靠性的理论基础。
定义1.1.1（适定性）的精确解读
定义原文：称初值问题(1.1.1)-(1.1.2)对初值\(y_a\)是适定的，如果存在常数\(K > 0\)和\(\eta > 0\)，使得对于任意\(0 < \varepsilon \leq \eta\)和
初值问题
存在解且满足\(|y(x) - z(x)| \leq K\varepsilon\)，其中\(y(x)\)是问题(1.1.1)-(1.1.2)的解。
深度解读：
- 适定性的本质：微分方程的解对输入数据（初始值\(y_a\)和右端函数\(f(x,y)\)）的连续依赖性。
- 物理意义：如果一个数学模型是适定的，那么输入数据的微小误差只会导致解的微小误差。这是所有科学计算的基本要求——否则，测量误差或计算误差会被无限放大，导致结果毫无意义。
- 定义中的两个扰动：
- 初始值的扰动：\(y_a \to \tilde{y}_a\)
- 右端函数的扰动：\(f(x,y) \to \tilde{f}(x,y)\)
这对应了实际问题中两种最常见的误差来源：初始条件的测量误差和模型本身的近似误差。
- 常数\(K\)的意义：称为条件数，它衡量了输入误差被放大的程度。\(K\)越小，问题的条件越好；\(K\)越大，问题的条件越差（病态）。
定理1.1.2（适定性定理）的完整证明
定理原文：如果方程(1.1.1)中的右端函数\(f = f(x,y)\)在区域\(\Omega\)上满足Lipschitz条件(1.1.3)，则初值问题(1.1.1)-(1.1.2)对任何初值\(y_a\)都是适定的。
证明：
设\(y(x)\)是原问题的解，\(z(x)\)是扰动问题的解。它们分别满足积分方程：
两式相减得：
将被积函数改写为：
两边取绝对值，并利用三角不等式和Lipschitz条件：
令\(C = \varepsilon [1 + (b - a)]\)，\(u(x) = |y(x) - z(x)|\)，则上式变为：
这是Gronwall不等式的一般形式。我们来求解这个不等式。
令\(v(x) = C + L \int_a^x u(t) dt\)，则\(v'(x) = L u(x)\)，且\(v(a) = C\)。由\(u(x) \leq v(x)\)得：
即
两边从\(a\)到\(x\)积分：
因此：
令\(K = [1 + (b - a)] e^{L(b - a)}\)，则：
这就证明了初值问题是适定的。
适定性的重要意义
- 数值方法的前提：只有当原问题是适定的，数值解法才有意义。如果原问题是不适定的，那么无论使用多么精确的数值方法，都无法得到可靠的结果。
- 误差分析的基础：适定性定理给出了输入误差与输出误差之间的定量关系，是数值方法误差分析的基础。
- 模型验证的标准：一个好的数学模型必须是适定的。如果一个模型是不适定的，那么它很可能没有正确反映物理现象的本质。
三、数值方法的基本思想：离散化
教材中指出："建立数值方法的过程也就是通过一些手段将问题(1.1.1)-(1.1.2)转化为在给定的有限个离散点\(\{x_m\}\)上近似(1.1.1)的有限差分或有限元方程的初值问题，这个过程通常称为离散化。"
离散化的具体步骤
- 区间离散化：将求解区间\([a,b]\)分成\(N\)个等距的子区间，分点为：\[x_0 = a, \quad x_1 = a + h, \quad x_2 = a + 2h, \quad \dots, \quad x_N = a + Nh = b \]其中\(h = (b - a)/N\)称为网格步长。
- 函数离散化：我们的目标不是求连续函数\(y = y(x)\)，而是求它在离散点\(x_0, x_1, \dots, x_N\)上的近似值\(y_0, y_1, \dots, y_N\)，其中\(y_0 = y_a\)是已知的初始值。
- 方程离散化：将微分方程\(y' = f(x,y)\)在每个离散点\(x_m\)上用某种近似方法转化为代数方程，从而建立从\(y_m\)计算\(y_{m+1}\)的递推公式。
三种主要的离散化方法
教材中提到了三种离散化方法，这里我简要介绍它们的本质：
| 离散化方法 | 基本思想 | 代表方法 |
|---|---|---|
| 直接化微商为差商法 | 用差商近似代替微商：\(y'(x_m) \approx \frac{y_{m+1} - y_m}{h}\) | 欧拉法 |
| Taylor级数展开法 | 将\(y(x_{m+1})\)在\(x_m\)处展开为Taylor级数，取前若干项作为近似 | 泰勒级数法、龙格-库塔法 |
| 数值积分方法 | 将微分方程从\(x_m\)到\(x_{m+1}\)积分，用数值积分公式近似计算积分 | 梯形法、亚当斯法 |
这三种方法本质上是等价的，它们都可以推导出相同的数值公式，但各自的出发点不同，适用于不同的场合。
四、数值方法的分类
教材将解常微分方程初值问题的数值方法分为两类：单步法和多步法。
1. 单步法
定义：计算\(y(x)\)在\(x = x_{m+1}\)处的值时仅用到\(x_m\)处的应变量及其导数值。
数学形式：
其中\(\phi(x, y, h)\)称为增量函数。
代表方法：
- 欧拉法（1.2节）
- 龙格-库塔法（1.3节）
优点：
- 自启动：只需要初始值\(y_0\)就可以开始计算
- 容易改变步长：在计算过程中可以根据需要随时调整步长
- 程序实现简单
缺点：
- 为了达到较高的精度，需要计算多个函数值（如四阶龙格-库塔法每步需要计算4个函数值）
- 对于某些问题，效率不如多步法
2. 多步法
定义：计算\(y(x)\)在\(x = x_{m+1}\)处的值时需要应变量及其导数在\(x_{m+1}\)左侧的多个网格结点处的值。
数学形式（k步法）：
其中\(\alpha_i\)和\(\beta_i\)是常数。如果\(\beta_{-1} = 0\)，则为显式多步法；如果\(\beta_{-1} \neq 0\)，则为隐式多步法。
代表方法：
- 线性多步方法（1.4节）
- 亚当斯显式法（Adams-Bashforth）
- 亚当斯隐式法（Adams-Moulton）
优点：
- 每步只需要计算一个新的函数值，效率高
- 容易构造高精度的方法
缺点：
- 不是自启动的：需要用单步法计算前k个初始值
- 改变步长比较麻烦
五、核心知识点归纳总结
我将1.1节所有内容（包括上一讲和本讲）归纳为以下表格，便于大家系统掌握：
| 知识点类别 | 具体内容 | 核心要点与公式 |
|---|---|---|
| 基本问题 | 一阶常微分方程初值问题 | \(\frac{dy}{dx} = f(x,y), \quad x \in [a,b]\) \(y(a) = y_a\) |
| 理论基础1 | 存在唯一性定理（全局版本） | 条件：\(f\)在\(\Omega = [a,b] \times \mathbb{R}\)上连续且关于\(y\)满足Lipschitz条件 结论：在\([a,b]\)上存在唯一解\(y(x) \in C^1[a,b]\) 证明方法：皮卡迭代法 |
| 理论基础2 | 适定性（连续依赖性） | 定义：输入的微小误差导致解的微小误差 定理：Lipschitz条件蕴含适定性 误差估计：\(|y(x) - z(x)| \leq K\varepsilon\)，其中\(K = [1 + (b-a)]e^{L(b-a)}\) |
| 数值方法基本思想 | 离散化 | 1. 区间离散化：\(x_m = a + mh, \quad h = (b-a)/N\) 2. 函数离散化：求\(y_m \approx y(x_m)\) 3. 方程离散化：建立递推公式 |
| 离散化方法 | 三种主要方法 | 1. 差商法：\(y'(x_m) \approx \frac{y_{m+1} - y_m}{h}\) 2. Taylor级数展开法 3. 数值积分法 |
| 数值方法分类 | 单步法 | 定义：仅用到前一个点的信息 形式：\(y_{m+1} = y_m + h\phi(x_m, y_m, h)\) 代表：欧拉法、龙格-库塔法 优点：自启动、易变步长 |
| 多步法 | 定义：用到前多个点的信息 形式：\(y_{m+1} = \sum_{i=0}^{k-1} \alpha_i y_{m-i} + h \sum_{i=-1}^{k-1} \beta_i f_{m-i}\) 代表：亚当斯法 优点：效率高、易构造高精度方法 |
|
| 关键概念 | Lipschitz条件 | \(|f(x,y) - f(x,z)| \leq L|y-z|\) 验证方法：若\(f_y\)有界，则\(L = \max |f_y|\) |
| Gronwall不等式 | 若\(u(x) \leq C + L \int_a^x u(t) dt\)，则\(u(x) \leq C e^{L(x-a)}\) |
六、学习建议
- 重点掌握两个核心定理：存在唯一性定理和适定性定理是整个常微分方程数值解法的理论基石，必须理解其条件、结论和证明过程。
- 深刻理解离散化思想：数值解法的本质是将连续问题转化为离散问题，这是贯穿整个数值分析学科的核心思想。
- 注意区分单步法和多步法：了解它们各自的优缺点和适用范围，为后续学习具体方法打下基础。
- 重视Gronwall不等式：它是微分方程理论中最重要的不等式之一，几乎所有的误差估计和稳定性分析都要用到它。
1.2 Euler方法 完整讲解与推导
作为拥有60多年教学经验的资深高级研究员，我将从三种推导视角、几何意义、误差分析、稳定性分析、改进方法和代码实现六个维度，对Euler方法进行系统、深入的讲解，补充教材中省略的所有关键证明和细节。
1.2.1 Euler方法及其几何意义
一、Euler方法的三种等价推导
Euler方法是求解常微分方程初值问题最简单、最基础的数值方法，它可以从三个完全不同的角度推导出来，这三个角度也代表了构造数值方法的三种基本思想。
1. 差商近似法（教材方法）
对于初值问题：
根据导数的定义，函数\(y(x)\)在点\(x_m\)处的导数为：
当步长\(h\)很小时，我们可以用向前差商近似代替导数：
代入微分方程\(y'(x_m) = f(x_m, y(x_m))\)，得：
整理得：
如果我们用\(y_m\)表示\(y(x_m)\)的近似值，就得到了显式Euler公式：
其中\(x_m = a + m h\)，\(m = 0, 1, 2, \dots, N-1\)，\(h = (b-a)/N\)为步长。
2. Taylor级数展开法
将\(y(x_{m+1})\)在点\(x_m\)处展开为Taylor级数：
其中\(\xi_m \in (x_m, x_{m+1})\)。
由于\(y'(x_m) = f(x_m, y(x_m))\)，代入得：
如果我们忽略高阶小项\(\frac{h^2}{2} y''(\xi_m)\)，就得到了Euler公式：
从这个推导可以清楚地看到，Euler方法的误差来源于Taylor级数的截断，这就是截断误差的由来。
3. 数值积分法
将微分方程\(y' = f(x,y)\)从\(x_m\)到\(x_{m+1}\)积分，得：
即：
如果我们用左矩形公式近似计算右边的积分：
代入(1.2.2)式，就得到了Euler公式。
三种推导方法的意义：
- 差商近似法直观易懂，体现了"以直代曲"的基本思想
- Taylor级数展开法揭示了截断误差的来源，是构造高阶数值方法的基础
- 数值积分法为构造更精确的数值方法（如梯形法、辛普森法）提供了统一的框架
二、Euler方法的几何意义
Euler方法有着非常直观的几何意义，因此也被称为折线法。
微分方程\(y' = f(x,y)\)的解是\((x,y)\)平面上的一族积分曲线，积分曲线上任意一点\((x,y)\)处的切线斜率等于\(f(x,y)\)。过点\((x_0, y_0)\)的积分曲线就是初值问题的解\(y = y(x)\)。
Euler方法的几何解释如下：
- 从初始点\(P_0(x_0, y_0)\)出发，作积分曲线在该点的切线，切线斜率为\(f(x_0, y_0)\)
- 该切线在\(x = x_1\)处与直线\(x = x_1\)交于点\(P_1(x_1, y_1)\)，其中\(y_1 = y_0 + h f(x_0, y_0)\)
- 再从点\(P_1(x_1, y_1)\)出发，作斜率为\(f(x_1, y_1)\)的切线，与直线\(x = x_2\)交于点\(P_2(x_2, y_2)\)
- 重复这个过程，得到一条折线\(P_0 P_1 P_2 \dots P_N\)
这条折线就是Euler方法得到的数值解，它近似代替了真实的积分曲线\(y = y(x)\)。
图1.2.1 Euler方法的几何说明
从图中可以清楚地看到：
- 每一步我们都用切线近似代替了积分曲线
- 随着步数的增加，折线会逐渐偏离真实的积分曲线，这就是误差积累的过程
- 步长\(h\)越小，切线与积分曲线的偏差就越小，数值解就越精确
1.2.2 Euler方法的误差分析
误差分析是数值方法的核心内容，它决定了数值方法的精度和可靠性。我们将Euler方法的误差分为局部截断误差和全局截断误差两类。
一、局部截断误差
定义：假设前一步的计算值是精确的，即\(y_m = y(x_m)\)，那么用数值方法计算\(y_{m+1}\)产生的误差称为局部截断误差，记为\(T_{m+1}\)：
Euler方法局部截断误差的推导：
由Taylor级数展开：
而Euler方法的计算值为：
两式相减得：
结论：
- Euler方法的局部截断误差为\(O(h^2)\)，即局部截断误差与步长\(h\)的平方成正比
- 我们称一个数值方法是p阶方法，如果它的局部截断误差为\(O(h^{p+1})\)
- 因此，Euler方法是一阶方法
二、全局截断误差
定义：在没有舍入误差的情况下，数值解\(y_m\)与精确解\(y(x_m)\)之间的误差称为全局截断误差，记为\(e_m\)：
全局截断误差是每一步局部截断误差积累的结果，它与局部截断误差有着密切的关系。
Euler方法全局截断误差的估计：
定理：如果函数\(f(x,y)\)关于\(y\)满足Lipschitz条件，Lipschitz常数为\(L\)，且存在常数\(M > 0\)，使得对所有\(x \in [a,b]\)，有\(|y''(x)| \leq M\)，则Euler方法的全局截断误差满足：
证明：
由局部截断误差的定义：
而Euler方法的计算值为：
两式相减得：
两边取绝对值，并利用Lipschitz条件和\(|T_{m+1}| \leq \frac{h^2 M}{2}\)：
这是一个递推不等式。我们用数学归纳法证明(1.2.3)式。
- 当\(m=0\)时，\(e_0 = y(x_0) - y_0 = 0\)，结论成立。
- 假设当\(m=k\)时结论成立，即：\[|e_k| \leq \frac{h M}{2 L} \left[ e^{L(x_k - a)} - 1 \right] \]那么当\(m=k+1\)时：\[\begin{align*} |e_{k+1}| &\leq (1 + h L) |e_k| + \frac{h^2 M}{2} \\ &\leq (1 + h L) \cdot \frac{h M}{2 L} \left[ e^{L(x_k - a)} - 1 \right] + \frac{h^2 M}{2} \\ &= \frac{h M}{2 L} \left[ (1 + h L) e^{L(x_k - a)} - (1 + h L) + h L \right] \\ &= \frac{h M}{2 L} \left[ (1 + h L) e^{L(x_k - a)} - 1 \right] \end{align*}\]由于\(1 + h L \leq e^{h L}\)（对所有\(h L \geq 0\)成立），因此：\[|e_{k+1}| \leq \frac{h M}{2 L} \left[ e^{h L} \cdot e^{L(x_k - a)} - 1 \right] = \frac{h M}{2 L} \left[ e^{L(x_{k+1} - a)} - 1 \right] \]结论成立。
由数学归纳法，(1.2.3)式对所有\(m = 0, 1, 2, \dots, N\)成立。
结论：
- Euler方法的全局截断误差为\(O(h)\)，即全局截断误差与步长\(h\)成正比
- 这与局部截断误差\(O(h^2)\)是一致的，因为全局误差是\(N = (b-a)/h\)个局部误差的积累，\(N \cdot O(h^2) = O(h)\)
- 当\(h \to 0\)时，全局截断误差趋向于0，因此Euler方法是收敛的
1.2.3 改进的Euler方法
Euler方法虽然简单，但精度较低（一阶方法）。为了提高精度，我们可以对Euler方法进行改进。
一、隐式Euler方法
如果我们用向后差商近似代替导数：
代入微分方程\(y'(x_{m+1}) = f(x_{m+1}, y(x_{m+1}))\)，得：
这就是隐式Euler公式（也称为向后Euler公式）。
特点：
- 公式(1.2.4)中\(y_{m+1}\)出现在等式两边，因此是隐式公式
- 一般情况下，需要通过迭代法求解\(y_{m+1}\)
- 局部截断误差：用Taylor级数展开可以证明，隐式Euler方法的局部截断误差也是\(O(h^2)\)，因此也是一阶方法
二、梯形方法
如果我们用梯形公式近似计算积分(1.2.2)：
代入(1.2.2)式，得到梯形公式：
局部截断误差分析：
将\(y(x_{m+1})\)和\(y'(x_{m+1})\)在\(x_m\)处展开为Taylor级数：
代入梯形公式的右边：
与\(y(x_{m+1})\)的展开式相减，得局部截断误差：
结论：
- 梯形方法的局部截断误差为\(O(h^3)\)，因此是二阶方法，精度比Euler方法高一阶
- 梯形公式也是隐式公式，需要迭代求解
三、预测-校正Euler方法（改进Euler方法）
隐式公式虽然精度较高，但需要迭代求解，计算量较大。为了兼顾精度和效率，我们通常采用预测-校正法：
- 预测步：用显式Euler公式计算一个初始近似值\(\bar{y}_{m+1}\)（预测值）\[\bar{y}_{m+1} = y_m + h f(x_m, y_m) \]
- 校正步：用预测值\(\bar{y}_{m+1}\)代入梯形公式的右边，计算校正值\(y_{m+1}\)\[y_{m+1} = y_m + \frac{h}{2} \left[ f(x_m, y_m) + f(x_{m+1}, \bar{y}_{m+1}) \right] \]
这就是预测-校正Euler方法，也称为改进Euler方法。
特点：
- 每步只需要计算两次函数值，计算量小
- 局部截断误差为\(O(h^3)\)，是二阶方法
- 不需要迭代，程序实现简单
1.2.4 Euler方法的稳定性分析
稳定性是数值方法的另一个重要性质，它研究的是计算过程中产生的舍入误差是否会被放大。
一、绝对稳定性的定义
我们用试验方程来分析数值方法的稳定性：
其中\(\lambda\)是复常数，\(\text{Re}(\lambda) < 0\)（保证微分方程本身是稳定的）。
定义：一个数值方法称为是绝对稳定的，如果用该方法求解试验方程(1.2.6)时，当步长\(h\)满足一定条件，数值解\(y_m\)满足：
使得数值方法绝对稳定的\(h\lambda\)的取值范围称为该方法的绝对稳定区域。
二、显式Euler方法的绝对稳定性
将显式Euler公式应用于试验方程(1.2.6)：
递推得：
要使\(\lim_{m \to \infty} y_m = 0\)，必须满足：
结论：
- 显式Euler方法的绝对稳定区域是复平面上以\((-1, 0)\)为圆心，半径为1的圆内部
- 对于实的\(\lambda < 0\)，绝对稳定区间为\(-2 < h \lambda < 0\)，即\(0 < h < -\frac{2}{\lambda}\)
三、隐式Euler方法的绝对稳定性
将隐式Euler公式应用于试验方程(1.2.6)：
整理得：
递推得：
要使\(\lim_{m \to \infty} y_m = 0\)，必须满足：
即：
结论：
- 隐式Euler方法的绝对稳定区域是复平面上以\((1, 0)\)为圆心，半径为1的圆外部
- 对于实的\(\lambda < 0\)，对任意\(h > 0\)都有\(|1 - h \lambda| = 1 - h \lambda > 1\)，因此隐式Euler方法对任意步长\(h\)都是绝对稳定的
- 这种对所有\(\text{Re}(\lambda) < 0\)都绝对稳定的方法称为A-稳定方法，特别适合求解刚性问题
1.2.5 数值例子与Python代码实现
一、数值例子
考虑初值问题：
该问题的精确解为\(y = \sqrt{1 + 2x}\)。
我们分别用显式Euler方法、隐式Euler方法和改进Euler方法求解，取步长\(h = 0.1\)，计算结果如下：
| \(x_m\) | 精确解\(y(x_m)\) | 显式Euler\(y_m\) | 显式Euler误差 | 改进Euler\(y_m\) | 改进Euler误差 |
|---|---|---|---|---|---|
| 0.0 | 1.0000 | 1.0000 | 0.0000 | 1.0000 | 0.0000 |
| 0.1 | 1.0954 | 1.1000 | 0.0046 | 1.0959 | 0.0005 |
| 0.2 | 1.1832 | 1.1918 | 0.0086 | 1.1841 | 0.0009 |
| 0.3 | 1.2649 | 1.2774 | 0.0125 | 1.2662 | 0.0013 |
| 0.4 | 1.3416 | 1.3582 | 0.0166 | 1.3434 | 0.0018 |
| 0.5 | 1.4142 | 1.4351 | 0.0209 | 1.4164 | 0.0022 |
| 0.6 | 1.4832 | 1.5090 | 0.0258 | 1.4860 | 0.0028 |
| 0.7 | 1.5492 | 1.5803 | 0.0311 | 1.5525 | 0.0033 |
| 0.8 | 1.6125 | 1.6498 | 0.0373 | 1.6165 | 0.0040 |
| 0.9 | 1.6733 | 1.7178 | 0.0445 | 1.6782 | 0.0049 |
| 1.0 | 1.7321 | 1.7848 | 0.0527 | 1.7379 | 0.0058 |
从计算结果可以清楚地看到：
- 改进Euler方法的精度比显式Euler方法高得多
- 随着\(x\)的增加，误差逐渐积累
- 显式Euler方法的误差约为改进Euler方法的10倍，这与它们的阶数一致（一阶vs二阶）
二、Python代码实现
import numpy as np
import matplotlib.pyplot as plt
def euler_explicit(f, a, b, y0, h):
"""
显式Euler方法求解常微分方程初值问题
参数:
f: 微分方程右端函数f(x,y)
a: 区间左端点
b: 区间右端点
y0: 初始值
h: 步长
返回:
x: 离散点数组
y: 数值解数组
"""
n = int((b - a) / h)
x = np.linspace(a, b, n+1)
y = np.zeros(n+1)
y[0] = y0
for m in range(n):
y[m+1] = y[m] + h * f(x[m], y[m])
return x, y
def euler_improved(f, a, b, y0, h):
"""
改进Euler方法（预测-校正法）求解常微分方程初值问题
参数:
f: 微分方程右端函数f(x,y)
a: 区间左端点
b: 区间右端点
y0: 初始值
h: 步长
返回:
x: 离散点数组
y: 数值解数组
"""
n = int((b - a) / h)
x = np.linspace(a, b, n+1)
y = np.zeros(n+1)
y[0] = y0
for m in range(n):
# 预测步
y_pred = y[m] + h * f(x[m], y[m])
# 校正步
y[m+1] = y[m] + h/2 * (f(x[m], y[m]) + f(x[m+1], y_pred))
return x, y
# 定义微分方程右端函数
def f(x, y):
return y - 2*x / y
# 精确解
def exact_solution(x):
return np.sqrt(1 + 2*x)
# 计算参数
a = 0.0
b = 1.0
y0 = 1.0
h = 0.1
# 计算数值解
x_euler, y_euler = euler_explicit(f, a, b, y0, h)
x_improved, y_improved = euler_improved(f, a, b, y0, h)
x_exact = np.linspace(a, b, 100)
y_exact = exact_solution(x_exact)
# 绘制结果
plt.figure(figsize=(10, 6))
plt.plot(x_exact, y_exact, 'k-', label='精确解')
plt.plot(x_euler, y_euler, 'ro--', label='显式Euler方法')
plt.plot(x_improved, y_improved, 'bs--', label='改进Euler方法')
plt.xlabel('x')
plt.ylabel('y')
plt.title('Euler方法求解常微分方程')
plt.legend()
plt.grid(True)
plt.show()
# 输出误差
print("x\t精确解\t\t显式Euler\t显式误差\t改进Euler\t改进误差")
print("-" * 70)
for i in range(len(x_euler)):
exact = exact_solution(x_euler[i])
err_euler = abs(y_euler[i] - exact)
err_improved = abs(y_improved[i] - exact)
print(f"{x_euler[i]:.1f}\t{exact:.4f}\t\t{y_euler[i]:.4f}\t\t{err_euler:.4f}\t\t{y_improved[i]:.4f}\t\t{err_improved:.4f}")
1.2.6 核心知识点归纳总结
| 方法名称 | 计算公式 | 局部截断误差 | 方法阶数 | 绝对稳定区域（实\(\lambda<0\)） | 显/隐式 | 每步函数计算次数 | 优点 | 缺点 |
|---|---|---|---|---|---|---|---|---|
| 显式Euler方法 | \(y_{m+1}=y_m+hf(x_m,y_m)\) | \(\frac{h^2}{2}y''(\xi_m)\) | 1阶 | \(0<h<-\frac{2}{\lambda}\) | 显式 | 1次 | 简单、自启动、易变步长 | 精度低、稳定性差 |
| 隐式Euler方法 | \(y_{m+1}=y_m+hf(x_{m+1},y_{m+1})\) | \(-\frac{h^2}{2}y''(\xi_m)\) | 1阶 | 任意\(h>0\)（A-稳定） | 隐式 | 迭代多次 | 稳定性好 | 精度低、需要迭代 |
| 梯形方法 | \(y_{m+1}=y_m+\frac{h}{2}[f(x_m,y_m)+f(x_{m+1},y_{m+1})]\) | \(-\frac{h^3}{12}y'''(\xi_m)\) | 2阶 | \(0<h<\infty\)（A-稳定） | 隐式 | 迭代多次 | 精度高、稳定性好 | 需要迭代 |
| 改进Euler方法 | 预测：\(\bar{y}_{m+1}=y_m+hf(x_m,y_m)\) 校正：\(y_{m+1}=y_m+\frac{h}{2}[f(x_m,y_m)+f(x_{m+1},\bar{y}_{m+1})]\) |
\(O(h^3)\) | 2阶 | \(0<h<-\frac{2}{\lambda}\) | 显式 | 2次 | 精度较高、不需要迭代、自启动 | 稳定性一般 |
学习建议
- 重点掌握三种推导方法：差商近似、Taylor展开和数值积分是构造所有数值方法的基本思想
- 深刻理解误差分析：局部截断误差决定了方法的阶数，全局截断误差决定了方法的收敛性
- 重视稳定性分析：稳定性是数值方法实用化的关键，隐式方法通常比显式方法更稳定
- 通过编程实践加深理解：亲手实现Euler方法和改进Euler方法，观察步长对误差的影响
1.2.2 Euler方法的误差分析 深度讲解与完整证明
作为拥有60多年教学经验的资深高级研究员，我将对教材中这部分内容进行逐行拆解式讲解，补充证明中每一步的理论依据，纠正常见理解误区，并完整推导教材中省略的全局截断误差估计。
一、局部截断误差的积分定义
教材从数值积分的角度重新定义了Euler方法的局部截断误差，这与我们之前从Taylor级数展开得到的定义是完全等价的，但更具一般性，是构造所有线性多步法的基础。
1. 积分形式的初值问题
对于初值问题
将方程在区间\([x_m, x_{m+1}]\)上积分，得到等价的积分形式：
这是一个精确等式，没有任何近似。所有基于数值积分的微分方程解法（包括Euler法、梯形法、亚当斯法等）都源于这个等式。
2. 左矩形积分与Euler方法
如果我们用左矩形求积公式近似计算(1.2.2)式右边的积分：
代入(1.2.2)式，得到：
如果用\(y_m\)代替\(y(x_m)\)，就得到了Euler公式：
3. 局部截断误差的精确定义
定义：假设前一步的计算值是精确的，即\(y_m = y(x_m)\)，那么用数值方法计算\(y_{m+1}\)产生的误差称为局部截断误差，记为\(R_m\)：
与Taylor展开定义的等价性：
将\(y(x_{m+1})\)在\(x_m\)处展开为Taylor级数：
与(1.2.4)式比较，立即得到：
这与我们之前得到的结果完全一致。
二、引理1.2.1的逐行证明与深度解读
教材给出的引理1.2.1是局部截断误差的一个先验界估计，它不依赖于解的二阶导数\(y''(x)\)，而是直接用右端函数\(f(x,y)\)的Lipschitz常数来表示，这在理论分析中非常有用。
引理1.2.1原文
设\(f(x,y)\)关于\(x\)和\(y\)均满足Lipschitz条件，\(K\)和\(L\)是相应的Lipschitz常数，即：
- 关于\(x\)的Lipschitz条件：\(|f(x_1,y) - f(x_2,y)| \leq K |x_1 - x_2|\)
- 关于\(y\)的Lipschitz条件：\(|f(x,y_1) - f(x,y_2)| \leq L |y_1 - y_2|\)
Euler方法的局部截断误差\(R_m\)满足：
其中
完整证明过程（逐行解读）
证明：由局部截断误差的定义(1.2.4)出发：
将\(h f(x_m, y(x_m))\)写成积分形式：
因此：
两边取绝对值，并利用积分的绝对值不等式：\(|\int_a^b g(s) ds| \leq \int_a^b |g(s)| ds\)，得：
关键技巧：加一项减一项
为了分别应用关于\(x\)和\(y\)的Lipschitz条件，我们在被积函数中插入一项\(f(x_m, y(s))\)：
代入(1)式，并利用三角不等式：\(|a + b| \leq |a| + |b|\)，得：
第一步：估计第一个积分
对第一个积分应用关于\(x\)的Lipschitz条件：
因此：
由于\(s \in [x_m, x_{m+1}]\)，所以\(|s - x_m| = s - x_m\)，积分得：
第二步：估计第二个积分
对第二个积分应用关于\(y\)的Lipschitz条件：
因此：
第三步：估计\(|y(s) - y(x_m)|\)
对函数\(y(x)\)在区间\([x_m, s]\)上应用拉格朗日微分中值定理：存在\(\xi \in (x_m, s)\)，使得
两边取绝对值：
其中\(M = \max_{a \leq x \leq b} |y'(x)|\)。代入(4)式，得：
第四步：合并结果
将(3)式和(5)式代入(2)式，得：
证毕。
重要说明
- 条件的意义：引理要求\(f(x,y)\)关于\(x\)也满足Lipschitz条件，这比我们之前用的只关于\(y\)满足Lipschitz条件更强。这是因为我们现在是从积分的角度进行估计，需要控制\(f\)随\(x\)的变化。
- 常数的意义：
- \(K\)：衡量\(f\)随\(x\)变化的快慢
- \(L\)：衡量\(f\)随\(y\)变化的快慢
- \(M\)：衡量解的导数的最大值，也就是解的变化率
- 误差阶数：从估计式可以看出，局部截断误差是\(O(h^2)\)，这与Taylor展开得到的结果一致。
三、全局截断误差的完整证明（教材省略部分）
教材中说"有了局部截断误差的界之后，可进一步研究各时间步的局部误差的累积，即估计整体误差\(\varepsilon_m = y(x_m) - y_m\)的界"，但没有给出具体证明。这里我给出完整、严谨的证明。
全局截断误差的定义
定义：在没有舍入误差的情况下，数值解\(y_m\)与精确解\(y(x_m)\)之间的误差称为全局截断误差，记为\(\varepsilon_m\)：
全局误差估计定理
定理：设\(f(x,y)\)关于\(x\)和\(y\)均满足Lipschitz条件，\(K\)和\(L\)是相应的Lipschitz常数，\(M = \max_{a \leq x \leq b} |f(x, y(x))|\)，则Euler方法的全局截断误差满足：
证明：
由局部截断误差的定义：
而Euler方法的计算值为：
(6)式减(7)式，得：
即：
两边取绝对值，并利用三角不等式和关于\(y\)的Lipschitz条件：
由引理1.2.1，\(|R_m| \leq \frac{h^2}{2} (K + L M) =: R\)，代入(9)式，得：
这是一个线性递推不等式。我们用数学归纳法证明(1.2.6)式。
-
基例：当\(m=0\)时，\(\varepsilon_0 = y(x_0) - y_0 = 0\)，结论成立。
-
归纳假设：假设当\(m=k\)时结论成立，即：
\[|\varepsilon_k| \leq \frac{R}{h L} \left[ e^{L(x_k - a)} - 1 \right] \]（注意：\(\frac{R}{h L} = \frac{h (K + L M)}{2 L}\)，与(1.2.6)式一致）
-
归纳步骤：当\(m=k+1\)时，由(10)式：
\[\begin{align*} |\varepsilon_{k+1}| &\leq (1 + h L) |\varepsilon_k| + R \\ &\leq (1 + h L) \cdot \frac{R}{h L} \left[ e^{L(x_k - a)} - 1 \right] + R \\ &= \frac{R}{h L} \left[ (1 + h L) e^{L(x_k - a)} - (1 + h L) + h L \right] \\ &= \frac{R}{h L} \left[ (1 + h L) e^{L(x_k - a)} - 1 \right] \end{align*}\]
利用不等式\(1 + t \leq e^t\)（对所有\(t \geq 0\)成立），令\(t = h L\)，得：
因此：
由数学归纳法，(1.2.6)式对所有\(m = 0, 1, 2, \dots, N\)成立。
结论
- 全局截断误差是\(O(h)\)，即当步长\(h\)减半时，全局误差也近似减半。
- 全局误差是局部误差的累积，由于共有\(N = (b-a)/h\)步，所以\(N \cdot O(h^2) = O(h)\)。
- 当\(h \to 0\)时，全局误差趋向于0，因此Euler方法是收敛的。
四、两种误差的对比与常见误区
局部截断误差 vs 全局截断误差
| 误差类型 | 定义 | 前提条件 | 阶数 | 物理意义 |
|---|---|---|---|---|
| 局部截断误差 | \(R_m = y(x_{m+1}) - [y(x_m) + h f(x_m, y(x_m))]\) | 前一步计算值精确：\(y_m = y(x_m)\) | \(O(h^2)\) | 单步计算产生的误差 |
| 全局截断误差 | \(\varepsilon_m = y(x_m) - y_m\) | 无舍入误差 | \(O(h)\) | 所有步局部误差累积的结果 |
常见误区
- 误区一："局部截断误差是\(O(h^2)\)，所以全局截断误差也是\(O(h^2)\)"
- 错误原因：忽略了误差的累积效应。每一步的局部误差都会传递到后续所有步，并且可能被放大。
- 误区二："步长越小，计算结果越精确"
- 错误原因：当步长太小时，舍入误差会变得显著。总误差 = 截断误差 + 舍入误差，存在一个最优步长使得总误差最小。
- 误区三："只要方法是收敛的，计算结果就一定可靠"
- 错误原因：收敛性是当\(h \to 0\)时的性质，而实际计算中\(h\)是有限的。还需要考虑方法的稳定性。
五、核心知识点归纳总结
| 知识点 | 核心内容 | 关键公式 |
|---|---|---|
| 积分形式初值问题 | 微分方程等价于积分方程 | \(y(x_{m+1}) = y(x_m) + \int_{x_m}^{x_{m+1}} f(s, y(s)) ds\) |
| 局部截断误差 | 单步计算产生的误差（前一步精确） | \(R_m = \int_{x_m}^{x_{m+1}} f(s, y(s)) ds - h f(x_m, y(x_m))\) |
| 引理1.2.1 | 局部截断误差的先验界 | \(|R_m| \leq \frac{h^2}{2} (K + L M)\) |
| 全局截断误差 | 所有步局部误差累积的结果 | \(\varepsilon_m = y(x_m) - y_m\) |
| 全局误差估计 | 全局截断误差的界 | \(|\varepsilon_m| \leq \frac{h (K + L M)}{2 L} \left[ e^{L(x_m - a)} - 1 \right]\) |
| 误差阶数 | 局部：\(O(h^2)\)，全局：\(O(h)\) | |
| 收敛性 | 当\(h \to 0\)时，数值解收敛于精确解 |
学习建议
- 重点掌握证明中的关键技巧：加一项减一项是分析中最常用的技巧之一，必须熟练掌握。
- 深刻理解两种误差的区别与联系：局部误差是单步误差，全局误差是累积误差，这是数值分析中最基本的概念。
- 注意定理的条件：所有误差估计都依赖于Lipschitz条件，这是保证解存在唯一和方法收敛的关键。
- 通过数值实验验证误差阶数：取不同的步长计算数值解，观察全局误差随步长的变化规律，验证\(O(h)\)的结论。
定理1.2.2（Euler方法收敛性定理）完整讲解与深度证明
作为拥有多年教学经验的资深高级研究员，我将对这个数值分析中最核心的收敛性定理进行逐行拆解式讲解，补充教材中省略的所有逻辑细节，揭示收敛性的本质，并澄清常见的理解误区。
一、定理完整陈述与条件解读
定理1.2.2原文
如果\(f(x,y)\)关于\(x\)和\(y\)均满足Lipschitz条件，\(K\)和\(L\)为相应的Lipschitz常数，且当\(h \to 0\)时，\(y_0 \to y(x_0)\)，则Euler方法(1.2.1)的解\(\{y_m\}\)一致收敛于初值问题(1.1.1)-(1.1.2)的解，且整体截断误差\(\varepsilon_m\)满足估计式
如果\(y_0 = y(x_0) = y(a)\)，即\(\varepsilon_0 = 0\)，则由上式得
这说明Euler方法的整体截断误差与\(h\)同阶。
条件深度解读
- Lipschitz条件：这是保证微分方程解存在唯一和数值方法收敛的核心条件
- 关于\(x\)的Lipschitz条件：\(|f(x_1,y) - f(x_2,y)| \leq K|x_1 - x_2|\)，控制\(f\)随\(x\)的变化率
- 关于\(y\)的Lipschitz条件：\(|f(x,y_1) - f(x,y_2)| \leq L|y_1 - y_2|\)，控制\(f\)随\(y\)的变化率
- 初始条件收敛：\(y_0 \to y(x_0)\)当\(h \to 0\)，这意味着初始误差随着步长减小而趋向于0
- 一致收敛：对于所有\(x_m \in [a,b]\)，都有\(\lim_{h \to 0} y_m = y(x_m)\)，即收敛性在整个区间上是均匀的
结论深度解读
- 误差分解：估计式(1.2.6)将整体误差分解为两个部分
- 第一部分：\(e^{L(b-a)}|\varepsilon_0|\)，是初始误差的传播项
- 第二部分：\(\frac{h}{2}\left(M + \frac{K}{L}\right)\left(e^{L(b-a)} - 1\right)\)，是局部截断误差的累积项
- 误差阶数：当初始误差\(\varepsilon_0 = 0\)时，整体截断误差是\(O(h)\)，即Euler方法是一阶收敛方法
- 收敛性：当\(h \to 0\)时，整体误差趋向于0，因此数值解收敛于精确解
二、定理完整证明（逐行解读）
第一步：建立误差递推关系
由局部截断误差的定义，精确解满足：
而Euler方法的数值解满足：
将两式相减，得到误差递推式：
即：
第二步：对误差递推式取绝对值并放缩
对(1.2.7)式两边取绝对值，并利用三角不等式：
再利用\(f\)关于\(y\)的Lipschitz条件：
代入上式得：
第三步：展开递推不等式
我们将递推不等式(1.2.8)从\(k=1\)到\(k=m\)反复展开：
- 当\(k=1\)时：\(|\varepsilon_1| \leq (1 + h L)|\varepsilon_0| + |R_0|\)
- 当\(k=2\)时：\(|\varepsilon_2| \leq (1 + h L)|\varepsilon_1| + |R_1| \leq (1 + h L)^2|\varepsilon_0| + (1 + h L)|R_0| + |R_1|\)
- 当\(k=3\)时：\(|\varepsilon_3| \leq (1 + h L)|\varepsilon_2| + |R_2| \leq (1 + h L)^3|\varepsilon_0| + (1 + h L)^2|R_0| + (1 + h L)|R_1| + |R_2|\)
以此类推，对于一般的\(k=m\)，有：
第四步：利用局部截断误差的界
由引理1.2.1，所有局部截断误差都有统一的上界：
代入(1.2.9)式，得：
第五步：计算等比数列的和
(1.2.10)式中的求和项是一个等比数列，首项为1，公比为\((1 + h L)\)，共有\(m\)项。等比数列求和公式为：
令\(q = 1 + h L\)，则：
代入(1.2.10)式，得：
第六步：利用指数不等式进行放缩
我们使用一个非常重要的不等式：对所有\(t \geq 0\)，有\(1 + t \leq e^t\)。
证明：令\(\varphi(t) = e^t - t - 1\)，则\(\varphi'(t) = e^t - 1 \geq 0\)对所有\(t \geq 0\)成立，因此\(\varphi(t)\)在\([0, +\infty)\)上单调递增。又因为\(\varphi(0) = 0\)，所以\(\varphi(t) \geq 0\)对所有\(t \geq 0\)成立，即\(e^t \geq 1 + t\)。
令\(t = h L\)，则：
两边取\(m\)次方，得：
由于\(x_m = a + m h \leq b\)，所以\(m h \leq b - a\)，因此：
第七步：得到最终的误差估计式
将(1.2.12)式代入(1.2.11)式，得：
将\(R = \frac{h^2}{2}(K + L M)\)代入上式，化简：
因此：
这就是估计式(1.2.6)。
第八步：证明收敛性
当初始误差\(\varepsilon_0 = 0\)时，估计式变为：
右边是一个与\(h\)成正比的量，记为\(C h\)，其中\(C = \frac{1}{2}\left(M + \frac{K}{L}\right)\left(e^{L(b-a)} - 1\right)\)是与\(h\)无关的常数。因此：
当\(h \to 0\)时，\(|\varepsilon_m| \to 0\)，即数值解\(y_m\)收敛于精确解\(y(x_m)\)。
又因为这个估计式对所有\(x_m \in [a,b]\)都成立，所以收敛是一致收敛。
证毕。
三、定理的深度解读与重要推论
1. 整体误差与局部误差的关系
定理最核心的结论是：Euler方法的整体截断误差的阶比局部截断误差低一阶。
- 局部截断误差：\(O(h^2)\)
- 整体截断误差：\(O(h)\)
原因解释：
整体误差是\(m\)个局部误差的累积，而\(m = \frac{x_m - a}{h} = O\left(\frac{1}{h}\right)\)。因此：
这个结论具有普遍性：对于大多数数值方法，整体截断误差的阶都比局部截断误差低一阶。
2. 初始误差的影响
估计式(1.2.6)中的第一项\(e^{L(b-a)}|\varepsilon_0|\)是初始误差的传播项。
- 系数\(e^{L(b-a)}\)称为误差放大因子
- 如果\(L\)很大，误差放大因子会非常大，这意味着初始误差会被迅速放大
- 这就是刚性问题的本质：Lipschitz常数\(L\)很大，导致数值方法的稳定性变差
3. 收敛性的定义
定义：一个数值方法称为是收敛的，如果对于任意固定的\(x \in [a,b]\)，当\(h \to 0\)且\(x_m = a + m h \to x\)时，有\(y_m \to y(x)\)。
定理1.2.2表明：如果微分方程的右端函数\(f(x,y)\)满足Lipschitz条件，那么Euler方法是收敛的。
4. 收敛性与稳定性的区别
- 收敛性：研究当\(h \to 0\)时，数值解是否趋向于精确解，是一个极限性质
- 稳定性：研究当\(h\)固定时，误差在计算过程中是否会被放大，是一个有限步长性质
收敛的方法不一定稳定，稳定的方法也不一定收敛。一个好的数值方法必须同时满足收敛性和稳定性。
四、常见误区澄清
- 误区一："只要方法是收敛的，计算结果就一定精确"
- 错误原因：收敛性是当\(h \to 0\)时的性质，而实际计算中\(h\)是有限的。对于有限的\(h\)，误差可能很大。
- 误区二："步长越小，计算结果越精确"
- 错误原因：当步长太小时，舍入误差会变得显著。总误差 = 截断误差 + 舍入误差，存在一个最优步长使得总误差最小。
- 误区三："Lipschitz条件是收敛的必要条件"
- 错误原因：Lipschitz条件是收敛的充分条件，但不是必要条件。存在一些不满足Lipschitz条件的问题，Euler方法仍然收敛。
五、核心知识点归纳总结
| 知识点 | 核心内容 | 关键公式 |
|---|---|---|
| 定理1.2.2 | Euler方法的收敛性定理 | |
| 条件 | 1. \(f(x,y)\)关于\(x\)和\(y\)满足Lipschitz条件 2. 初始误差\(\varepsilon_0 \to 0\)当\(h \to 0\) |
|
| 结论 | 1. Euler方法一致收敛于精确解 2. 整体截断误差\(|\varepsilon_m| = O(h)\) |
|
| 误差估计式 | 整体误差分解为初始误差传播和局部误差累积 | \(|\varepsilon_m| \leq e^{L(b-a)}|\varepsilon_0| + \frac{h}{2}\left(M + \frac{K}{L}\right)\left(e^{L(b-a)} - 1\right)\) |
| 误差阶数 | 局部截断误差：\(O(h^2)\) 整体截断误差：\(O(h)\) |
|
| 关键不等式 | 指数不等式 | \(1 + t \leq e^t, \quad t \geq 0\) |
| 等比数列求和 | 递推不等式展开的关键 | \(\sum_{i=0}^{m-1} q^i = \frac{q^m - 1}{q - 1}\) |
| 收敛性定义 | 当\(h \to 0\)时，数值解趋向于精确解 | \(\lim_{h \to 0} y_m = y(x_m)\) |
学习建议
- 重点掌握递推不等式的解法：这是所有数值方法收敛性分析的基础，必须熟练掌握展开递推和等比数列求和的技巧。
- 深刻理解误差分解：整体误差由初始误差和局部误差两部分组成，这是理解数值方法行为的关键。
- 区分收敛性和稳定性：这是两个不同的概念，收敛性是极限性质，稳定性是有限步长性质。
- 通过数值实验验证：取不同的步长计算数值解，观察整体误差随步长的变化规律，验证\(O(h)\)的结论。
1.2.3 Euler方法的稳定性 教材内容逐行精讲与深度拓展
作为拥有60多年教学经验的资深高级研究员，我将严格按照教材逻辑，对这部分内容进行逐句拆解式讲解，补充证明中每一步的理论依据，澄清教材中隐含的概念，并拓展实际计算中至关重要的绝对稳定性内容。
一、定义1.2.1（零稳定性）的精确解读
教材原文
定义1.2.1 称Euler方法(1.2.1)是稳定的，如果存在常数\(C\)和\(h_0\)，使得Euler方法的解\(y_m\)和\(z_m\)满足
其中\(y_m\)和\(z_m\)是方程(1.2.1)分别以\(y_0\)和\(z_0\)为初始值的解。
该稳定性表示，对于任意\(h \in (0,h_0)\)，Euler方法(1.2.1)的精确解连续地依赖于初始值。
逐句深度解读
-
"存在常数\(C\)和\(h_0\)"
- \(C\)称为误差放大因子，它是一个与步长\(h\)和步数\(m\)无关的常数，只依赖于求解区间\([a,b]\)和右端函数\(f(x,y)\)的性质
- \(h_0\)称为稳定步长上限，只要步长\(h\)小于这个临界值，稳定性就成立
- 这是一个存在性条件，不要求我们找到最优的\(C\)和\(h_0\)，只要证明它们存在即可
-
"\(|y_m - z_m| \leq C|y_0 - z_0|\)"
- 左边是第\(m\)步的误差，右边是初始误差乘以一个常数
- 这意味着初始误差不会随着计算步数的增加而无限放大，最多被放大\(C\)倍
- 这是数值方法能够实用的最基本要求：如果初始的微小误差会被无限放大，那么计算结果将毫无意义
-
"\(0 < h < h_0, \ a \leq mh \leq b\)"
- 第一个条件说明稳定性是当步长趋向于0时的性质
- 第二个条件说明稳定性是在有限区间\([a,b]\)上成立的
- 这种稳定性在数值分析中被称为零稳定性（也叫李雅普诺夫稳定性）
-
"精确解连续地依赖于初始值"
- 这是微分方程本身适定性的离散版本
- 微分方程的适定性要求解连续依赖于初始值，而数值方法的稳定性要求数值解也连续依赖于初始值
- 只有这样，数值方法才能正确反映原微分方程的性质
二、定理1.2.3的完整证明与逐行推导
教材原文
定理1.2.3 在定理1.2.2的假设条件下，Euler方法(1.2.1)是稳定的。
证明 考虑
和
两式相减，并记\(e_m = y_m - z_m\)，则有
从而，当\(a \leq x_m = a + mh \leq b\)时，有\(|e_m| \leq e^{L(b-a)}|e_0| = C|e_0|\).
逐行证明与理论依据
步骤1：构造两个不同初始值的解
我们考虑Euler方法的两个解：
- \(y_m\)：以\(y_0\)为初始值的解
- \(z_m\)：以\(z_0\)为初始值的解
它们分别满足递推公式：
步骤2：定义误差并建立误差递推关系
将(1)式和(2)式相减，得到：
定义误差\(e_m = y_m - z_m\)，则上式变为：
步骤3：对误差递推式取绝对值并放缩
对(3)式两边取绝对值，并利用三角不等式\(|a + b| \leq |a| + |b|\)，得：
步骤4：应用Lipschitz条件
根据定理1.2.2的假设，\(f(x,y)\)关于\(y\)满足Lipschitz条件，即存在常数\(L > 0\)，使得：
将其代入(4)式，得到：
步骤5：展开递推不等式
(5)式是一个简单的递推不等式，我们可以将其逐次展开：
- 当\(m=0\)时：\(|e_1| \leq (1 + h L) |e_0|\)
- 当\(m=1\)时：\(|e_2| \leq (1 + h L) |e_1| \leq (1 + h L)^2 |e_0|\)
- 当\(m=2\)时：\(|e_3| \leq (1 + h L) |e_2| \leq (1 + h L)^3 |e_0|\)
- ...
- 一般地，对于任意\(m \geq 0\)，有：\[|e_m| \leq (1 + h L)^m |e_0| \tag{6} \]
步骤6：利用指数不等式进行最终放缩
我们使用数值分析中最基本的不等式：对所有\(t \geq 0\)，有\(1 + t \leq e^t\)。
证明：令\(\varphi(t) = e^t - t - 1\)，则\(\varphi'(t) = e^t - 1 \geq 0\)对所有\(t \geq 0\)成立，因此\(\varphi(t)\)在\([0, +\infty)\)上单调递增。又因为\(\varphi(0) = 0\)，所以\(\varphi(t) \geq 0\)对所有\(t \geq 0\)成立，即\(e^t \geq 1 + t\)。
令\(t = h L\)，则：
两边取\(m\)次方，得：
由于\(x_m = a + m h \leq b\)，所以\(m h \leq b - a\)，因此：
将(7)式和(8)式代入(6)式，得到：
令\(C = e^{L(b - a)}\)，则：
这就证明了Euler方法是稳定的。
证毕。
证明的关键要点总结
- 误差递推：通过两个解相减得到误差的递推关系，这是所有稳定性分析的标准方法
- Lipschitz条件：这是证明稳定性的核心条件，它控制了右端函数对误差的放大作用
- 指数不等式：\(1 + t \leq e^t\)是数值分析中使用最频繁的不等式，用于将乘积形式转化为指数形式
- 有限区间：稳定性在有限区间\([a,b]\)上成立，这是因为\(m h \leq b - a\)，所以指数项是有界的
三、教材隐含概念的澄清：零稳定性 vs 绝对稳定性
教材中只介绍了零稳定性，但在实际数值计算中，绝对稳定性更为重要。这两种稳定性描述的是完全不同的情况，必须严格区分。
1. 零稳定性（教材定义）
- 研究场景：当步长\(h \to 0\)时的误差传播
- 核心问题：当步长无限减小时，初始误差是否会被无限放大
- 理论意义：是证明数值方法收敛性的必要条件（Lax等价定理）
- Euler方法的情况：在Lipschitz条件下，Euler方法总是零稳定的
2. 绝对稳定性（实际计算必需）
- 研究场景：当步长\(h\)固定为有限值时的误差传播
- 核心问题：对于给定的有限步长\(h\)，计算过程中产生的舍入误差是否会被无限放大
- 实际意义：直接决定了数值方法能否在实际计算中使用
- Euler方法的情况：显式Euler方法只有当步长\(h\)小于某个临界值时才是绝对稳定的
3. 为什么绝对稳定性更重要？
在实际计算中，我们不可能取\(h \to 0\)，只能取有限的步长\(h\)。如果一个方法不是绝对稳定的，那么即使它是零稳定的，计算过程中产生的舍入误差也会被迅速放大，导致计算结果完全不可靠。
4. 显式Euler方法的绝对稳定性分析
我们用试验方程\(y' = \lambda y\)（\(\text{Re}(\lambda) < 0\)，保证微分方程本身稳定）来分析绝对稳定性。
将显式Euler公式应用于试验方程：
递推得：
要使\(\lim_{m \to \infty} y_m = 0\)（误差不被放大），必须满足：
结论：
- 显式Euler方法的绝对稳定区域是复平面上以\((-1, 0)\)为圆心，半径为1的圆内部
- 对于实的\(\lambda < 0\)，绝对稳定区间为：\[0 < h < -\frac{2}{\lambda} \]
这意味着，当\(\lambda\)的绝对值很大时（刚性问题），显式Euler方法必须使用非常小的步长才能保证稳定，这会导致计算量急剧增加。
四、稳定性与收敛性的关系：Lax等价定理
数值分析中最重要的定理之一是Lax等价定理，它揭示了收敛性与稳定性之间的本质联系。
Lax等价定理
对于一个相容的数值方法，收敛性等价于零稳定性。
概念解释
- 相容性：数值方法的局部截断误差趋向于0当\(h \to 0\)。对于Euler方法，局部截断误差是\(O(h^2)\)，因此是相容的。
- 零稳定性：就是教材中定义的稳定性。
- 收敛性：当\(h \to 0\)时，数值解趋向于精确解。
对Euler方法的应用
- 相容性：满足（局部截断误差\(O(h^2)\)）
- 零稳定性：满足（定理1.2.3）
- 因此，Euler方法是收敛的（定理1.2.2）
这就解释了为什么我们要先证明稳定性，再证明收敛性——对于相容的方法，稳定性是收敛性的充要条件。
五、核心知识点归纳总结
| 知识点 | 核心内容 | 关键公式/结论 |
|---|---|---|
| 零稳定性定义 | 初始误差不会随步数增加而无限放大 | \(|y_m - z_m| \leq C|y_0 - z_0|, \ 0 < h < h_0\) |
| 定理1.2.3 | 在Lipschitz条件下，Euler方法是零稳定的 | \(|e_m| \leq e^{L(b-a)}|e_0|\) |
| 证明关键步骤 | 1. 构造两个解并相减得到误差递推 2. 应用Lipschitz条件放缩 3. 展开递推不等式 4. 利用指数不等式得到最终界 |
\(|e_{m+1}| \leq (1 + hL)|e_m|\) \((1 + hL)^m \leq e^{L(b-a)}\) |
| 两种稳定性对比 | 零稳定性：\(h \to 0\)时的理论性质 绝对稳定性：\(h\)固定时的实际性质 |
显式Euler绝对稳定区间：\(0 < h < -\frac{2}{\lambda}\) |
| Lax等价定理 | 相容的数值方法，收敛性等价于零稳定性 |
学习建议
- 熟练掌握稳定性证明的标准流程：构造两个解→相减得到误差递推→应用Lipschitz条件→展开递推→指数不等式放缩，这是所有线性多步法稳定性证明的通用方法。
- 严格区分零稳定性和绝对稳定性：零稳定性是理论基础，绝对稳定性是实际计算的关键，不要混淆这两个概念。
- 理解Lax等价定理的深刻意义：它将数值方法的三个核心性质（相容性、稳定性、收敛性）联系在了一起，是数值分析的基石。
- 通过数值实验验证绝对稳定性：取不同的步长计算试验方程\(y' = -100y\)，观察显式Euler方法在步长超过临界值时的不稳定现象。
1.2.4 改进的Euler方法 完整讲解与深度证明
作为拥有多年教学经验的资深高级研究员，我将对这部分内容进行系统讲解，完整推导梯形公式的局部截断误差，严格证明迭代求解的收敛性，并补充实际计算中最常用的预测-校正Euler方法。
一、为什么需要改进Euler方法？
显式Euler方法虽然简单易实现，但存在一个致命的缺点：精度太低。
精度问题的直观体现
Euler方法的整体截断误差是\(O(h)\)，这意味着要将误差减小一个数量级，步长\(h\)必须减小一个数量级，而计算步数会增加一个数量级。
教材中给出了一个具体的例子：
- 要求误差达到\(O(10^{-6})\)
- 则需要\(h \approx 10^{-6}\)
- 在区间\([0,1]\)上需要计算\(10^6\)步
这在实际计算中是难以接受的，特别是对于大规模问题。因此，我们必须寻求精度更高的数值方法。
提高精度的基本思路
回顾Euler方法的推导过程：它是用左矩形求积公式近似计算积分
得到的。左矩形公式的精度是\(O(h^2)\)，因此Euler方法的局部截断误差是\(O(h^2)\)，整体截断误差是\(O(h)\)。
一个自然的想法是：用更高精度的数值积分公式代替左矩形公式，就可以得到更高精度的微分方程数值解法。
梯形求积公式是比左矩形公式高一阶的数值积分公式，它的精度是\(O(h^3)\)，因此用它来近似计算积分，就可以得到二阶精度的数值方法。
二、梯形公式的推导
积分形式的初值问题
对于初值问题
将方程在区间\([x_m, x_{m+1}]\)上积分，得到精确等式：
梯形求积公式
梯形求积公式是用梯形的面积近似代替曲边梯形的面积：
其中\(h = x_{m+1} - x_m\)是步长。
将梯形求积公式应用于(1.2.2)式右边的积分，得到：
如果用\(y_m\)代替\(y(x_m)\)，\(y_{m+1}\)代替\(y(x_{m+1})\)，就得到了梯形公式：
三、梯形公式的局部截断误差（完整证明）
教材原文
其中
容易验证，当\(y(x)\)足够光滑时，\(R_m^{(1)}\)的阶比\(R_m\)高，且有
完整证明过程
证明：由于\(y(x)\)是微分方程的解，所以\(f(x, y(x)) = y'(x)\)。因此，局部截断误差可以写为：
我们通过变量替换将积分区间变换到\([0,1]\)。令
则\(dx = h d\tau\)，当\(x = x_m\)时\(\tau = 0\)，当\(x = x_{m+1}\)时\(\tau = 1\)。代入(9)式，得：
现在，将\(y'(x_m + \tau h)\)在点\(x_m\)处展开为带积分余项的泰勒级数：
但教材中使用了一种更巧妙的展开方式：将\(y'(x)\)在区间\([x_m, x_{m+1}]\)上表示为线性插值加上余项。线性插值多项式为：
它在\(\tau = 0\)时等于\(y'(x_m)\)，在\(\tau = 1\)时等于\(y'(x_{m+1})\)。
根据插值余项定理，有：
将(11)式代入(10)式，得：
我们分别计算三个积分：
- 第一个积分：\[\int_0^1 y'(x_m) d\tau = y'(x_m) \]
- 第二个积分：\[\int_0^1 \tau \left(y'(x_{m+1}) - y'(x_m)\right) d\tau = \left(y'(x_{m+1}) - y'(x_m)\right) \cdot \left. \frac{\tau^2}{2} \right|_0^1 = \frac{1}{2} \left(y'(x_{m+1}) - y'(x_m)\right) \]
- 第三个积分：\[\int_0^1 \frac{h^2}{2} \tau (\tau - 1) y'''(x_m + \theta h) d\tau = \frac{h^2}{2} y'''(x_m + \xi h) \int_0^1 \tau (\tau - 1) d\tau \]其中\(\xi \in [0,1]\)（积分中值定理）。计算积分：\[\int_0^1 \tau (\tau - 1) d\tau = \int_0^1 (\tau^2 - \tau) d\tau = \left. \left( \frac{\tau^3}{3} - \frac{\tau^2}{2} \right) \right|_0^1 = \frac{1}{3} - \frac{1}{2} = -\frac{1}{6} \]
将三个积分的结果代入，得：
证毕。
结论
- 梯形公式的局部截断误差是\(O(h^3)\)，比Euler方法的\(O(h^2)\)高一阶
- 因此，梯形公式是二阶方法，整体截断误差是\(O(h^2)\)
- 这意味着要达到相同的精度，梯形方法可以使用比Euler方法大得多的步长，从而大大减少计算量
四、梯形公式的隐式特性与迭代求解
隐式公式的问题
梯形公式(1.2.12)的右边出现了未知量\(y_{m+1}\)，因此它是一个隐式公式。
- 当\(f(x,y)\)是\(y\)的线性函数时，可以直接解出\(y_{m+1}\)
- 当\(f(x,y)\)是\(y\)的非线性函数时，(1.2.12)是一个关于\(y_{m+1}\)的非线性方程，需要通过迭代法求解
迭代求解公式
我们采用不动点迭代法求解梯形公式。迭代公式为：
迭代收敛性证明
定理：如果\(f(x,y)\)关于\(y\)满足Lipschitz条件，Lipschitz常数为\(L\)，且步长\(h\)满足
则迭代公式(1.2.13)收敛。
证明：定义迭代函数
则迭代公式可以写为
根据压缩映射原理，如果\(G(y)\)是压缩映射，即存在常数\(0 \leq k < 1\)，使得
对所有\(y, z\)成立，则迭代收敛。
现在计算\(|G(y) - G(z)|\)：
当\(h L \leq 2\)时，\(\frac{h L}{2} \leq 1\)。特别地，当\(h L < 2\)时，\(\frac{h L}{2} < 1\)，因此\(G(y)\)是压缩映射，迭代收敛。
当\(h L = 2\)时，\(\frac{h L}{2} = 1\)，此时迭代仍然收敛（在实际计算中，通常取\(h L < 2\)以保证收敛速度）。
证毕。
迭代次数
在实际计算中，梯形公式的迭代收敛速度非常快，通常只需要2~3次迭代就可以达到足够的精度。这是因为初始猜测\(y_{m+1}^{[0]}\)已经是\(O(h^2)\)精度的近似值，而每次迭代可以将误差提高一个数量级。
五、预测-校正Euler方法（改进Euler方法）
梯形公式虽然精度高，但需要迭代求解，计算量较大。为了兼顾精度和效率，我们通常采用预测-校正法。
基本思想
- 预测步：用显式Euler公式计算一个初始近似值\(\bar{y}_{m+1}\)（预测值）
- 校正步：用预测值\(\bar{y}_{m+1}\)代替梯形公式右边的\(y_{m+1}\)，计算校正值\(y_{m+1}\)
公式
局部截断误差分析
可以证明，预测-校正Euler方法的局部截断误差也是\(O(h^3)\)，因此它也是二阶方法，精度与梯形公式相同。
证明思路：将预测值\(\bar{y}_{m+1}\)展开为泰勒级数，然后代入校正步公式，与\(y(x_{m+1})\)的泰勒级数展开式比较，可以发现它们的前三项完全相同，误差项是\(O(h^3)\)。
优点
- 是显式公式，不需要迭代，计算量小（每步只需要计算两次函数值）
- 是自启动的，只需要初始值\(y_0\)就可以开始计算
- 精度高（二阶方法），比显式Euler方法精度高一阶
缺点
- 绝对稳定区间与显式Euler方法相同，都是\(0 < h < -\frac{2}{\lambda}\)，是条件稳定的
- 对于刚性问题，仍然需要使用很小的步长
六、四种Euler类方法的全面对比
| 方法名称 | 计算公式 | 局部截断误差 | 方法阶数 | 显/隐式 | 每步函数计算次数 | 绝对稳定区间（实\(\lambda<0\)） | 核心优势 | 核心缺点 |
|---|---|---|---|---|---|---|---|---|
| 显式Euler | \(y_{m+1}=y_m+hf(x_m,y_m)\) | \(\frac{h^2}{2}y''(\xi_m)\) | 1阶 | 显式 | 1次 | \(0<h<-\frac{2}{\lambda}\) | 最简单、自启动 | 精度低、稳定性差 |
| 隐式Euler | \(y_{m+1}=y_m+hf(x_{m+1},y_{m+1})\) | \(-\frac{h^2}{2}y''(\xi_m)\) | 1阶 | 隐式 | 迭代多次 | 任意\(h>0\)（A-稳定） | 稳定性极好 | 精度低、需要迭代 |
| 梯形方法 | \(y_{m+1}=y_m+\frac{h}{2}[f(x_m,y_m)+f(x_{m+1},y_{m+1})]\) | \(-\frac{h^3}{12}y'''(\xi_m)\) | 2阶 | 隐式 | 迭代多次 | 任意\(h>0\)（A-稳定） | 精度高、稳定性极好 | 需要迭代 |
| 改进Euler | 预测：\(\bar{y}_{m+1}=y_m+hf(x_m,y_m)\) 校正：\(y_{m+1}=y_m+\frac{h}{2}[f(x_m,y_m)+f(x_{m+1},\bar{y}_{m+1})]\) |
\(O(h^3)\) | 2阶 | 显式 | 2次 | \(0<h<-\frac{2}{\lambda}\) | 精度较高、不需要迭代、自启动 | 稳定性一般 |
七、梯形公式的收敛性与稳定性
教材中说"可模仿定理1.2.2和定理1.2.3的证明建立梯形公式的收敛性和稳定性"，这里给出简要结论：
收敛性
在定理1.2.2的假设条件下，梯形公式是收敛的，且整体截断误差是\(O(h^2)\)。
稳定性
- 零稳定性：梯形公式是零稳定的
- 绝对稳定性：梯形公式是A-稳定的，对任意步长\(h > 0\)都是绝对稳定的
这是梯形公式最突出的优点之一，它既具有二阶精度，又具有极好的稳定性，是求解刚性问题的常用方法之一。
八、核心知识点归纳总结
| 知识点 | 核心内容 | 关键公式/结论 |
|---|---|---|
| 改进的动机 | Euler方法精度低（一阶），需要高阶方法 | 误差\(O(h)\)，要达到\(10^{-6}\)精度需要\(10^6\)步 |
| 梯形公式推导 | 用梯形求积公式近似积分 | \(y_{m+1}=y_m+\frac{h}{2}[f(x_m,y_m)+f(x_{m+1},y_{m+1})]\) |
| 局部截断误差 | 梯形公式是二阶方法 | \(R_m^{(1)}=-\frac{h^3}{12}y'''(x_m+\xi h)\) |
| 迭代求解 | 隐式公式需要迭代求解 | \(y_{m+1}^{[\nu+1]}=y_m+\frac{h}{2}[f(x_m,y_m)+f(x_{m+1},y_{m+1}^{[\nu]})]\) |
| 迭代收敛条件 | \(hL \leq 2\) | \(|G(y)-G(z)| \leq \frac{hL}{2}|y-z|\) |
| 预测-校正法 | 显式二阶方法，实际最常用 | 预测步用显式Euler，校正步用梯形公式 |
| 稳定性 | 梯形公式是A-稳定的 | 对任意\(h>0\)都绝对稳定 |
学习建议
- 重点掌握梯形公式的局部截断误差推导：这是理解高阶数值方法精度来源的关键。
- 理解隐式公式的迭代求解：隐式方法虽然需要迭代，但具有极好的稳定性，是求解刚性问题的基础。
- 熟练掌握预测-校正Euler方法：这是实际计算中最常用的二阶显式方法，兼顾了精度和效率。
- 通过数值实验对比不同方法：取相同的步长，计算同一个初值问题，观察不同方法的误差和稳定性差异。
1.3 Runge-Kutta方法 完整讲解与深度推导
作为拥有多年教学经验的资深高级研究员，我将对这部分内容进行系统、深入的讲解，完整推导阶数条件，揭示Runge-Kutta方法的本质，并详细介绍工程中最常用的各种公式。
一、Runge-Kutta方法的基本思想
1. 直接Taylor展开法的局限性
理论上，我们可以用Taylor级数展开法构造任意高阶的数值方法。对于初值问题
将\(y(x_{m+1})\)在\(x_m\)处展开为Taylor级数：
其中各阶导数可以通过微分方程逐次求导得到：
- 一阶导数：\(y'(x_m) = f(x_m, y(x_m))\)
- 二阶导数：\(y''(x_m) = f_x(x_m, y(x_m)) + f_y(x_m, y(x_m)) y'(x_m)\)
- 三阶导数：\(y'''(x_m) = f_{xx} + 2f_{xy} y' + f_{yy} (y')^2 + f_y y''\)
- ...
局限性：
- 高阶导数的计算非常繁琐，特别是对于复杂的\(f(x,y)\)
- 得到的公式形式复杂，不适合实际应用
- 每提高一阶精度，计算量会急剧增加
2. Runge的天才想法
1895年，德国数学家Carl Runge提出了一个革命性的思想：间接使用Taylor展开。
核心思想：
- 不直接计算\(y\)的高阶导数
- 而是用函数\(f\)在\(s\)个不同点处的值的线性组合来代替高阶导数
- 然后通过Taylor展开比较系数，确定这些线性组合的系数
优势：
- 避免了高阶导数的计算
- 公式形式简单，易于编程实现
- 可以构造出高精度的单步方法
二、显式Runge-Kutta方法的一般形式
定义1.3.1 完整解读
s级显式Runge-Kutta方法的一般形式为：
其中\(k_i\)满足：
参数含义：
- \(s\)：级数，表示每步需要计算函数\(f\)的次数
- \(c_i\)：节点参数，表示第\(i\)个点相对于\(x_m\)的横坐标偏移
- \(a_{i,j}\)：内部权重，表示计算\(k_i\)时\(k_j\)的权重
- \(b_i\)：外部权重，表示最终组合时\(k_i\)的权重
Butcher表（Butcher阵列）
1964年，Butcher提出了一种非常直观简洁的方式来表示Runge-Kutta方法，称为Butcher表：
解读：
- 第一列是节点参数\(c_i\)（\(c_1=0\)）
- 第\(i\)行（\(i \geq 2\)）是内部权重\(a_{i,1}, \dots, a_{i,i-1}\)
- 最后一行是外部权重\(b_i\)
- 由于是显式方法，矩阵\(A = (a_{i,j})\)是严格下三角矩阵
三、阶数条件的完整推导（以s=2为例）
1. 基本思路
要使s级Runge-Kutta方法达到p阶精度，必须满足：
具体步骤：
- 将\(y(x_{m+1})\)在\(x_m\)处展开为Taylor级数
- 将各个\(k_i\)在\(x_m\)处展开为Taylor级数
- 将\(k_i\)代入\(y_{m+1}\)的表达式，得到\(y_{m+1}\)的Taylor展开式
- 比较两个展开式中\(h^i\)（\(i=1,2,\dots,p\)）的系数，令其相等
- 得到确定系数\(a_{i,j}, b_i, c_i\)的代数方程组
2. 二级二阶Runge-Kutta方法的阶数条件推导
对于\(s=2\)的情况，方法形式为：
第一步：展开\(k_2\)
将二元函数\(f(x + \Delta x, y + \Delta y)\)在\((x_m, y_m)\)处展开为Taylor级数：
令\(\Delta x = c_2 h\)，\(\Delta y = h a_{2,1} k_1 = h a_{2,1} f_m\)，代入得：
第二步：展开\(y_{m+1}\)
将\(k_1\)和\(k_2\)代入\(y_{m+1}\)的表达式：
第三步：展开精确解\(y(x_{m+1})\)
第四步：比较系数
比较\(y_{m+1}\)和\(y(x_{m+1})\)的展开式中\(h\)和\(h^2\)项的系数，令其相等：
- \(h\)项系数：\(b_1 + b_2 = 1\)
- \(h^2\)项中\(f_{x,m}\)的系数：\(b_2 c_2 = \frac{1}{2}\)
- \(h^2\)项中\(f_m f_{y,m}\)的系数：\(b_2 a_{2,1} = \frac{1}{2}\)
因此，我们得到了二级二阶Runge-Kutta方法的阶数条件：
3. 解的自由度与典型方法
方程组(1.3.6)有3个方程，但有4个未知数\(b_1, b_2, c_2, a_{2,1}\)，因此有一个自由参数。我们可以取\(c_2\)为自由参数，然后解出其他三个参数：
典型的二级二阶方法：
-
中点公式（取\(c_2 = \frac{1}{2}\)）
\[b_1 = 0, \quad b_2 = 1, \quad a_{2,1} = \frac{1}{2} \]公式：
\[y_{m+1} = y_m + h f\left(x_m + \frac{h}{2}, y_m + \frac{h}{2} f_m\right) \]Butcher表：
\[\begin{array}{c|cc} 0 & & \\ \frac{1}{2} & \frac{1}{2} & \\ \hline & 0 & 1 \\ \end{array} \] -
Heun公式（取\(c_2 = \frac{2}{3}\)）
\[b_1 = \frac{1}{4}, \quad b_2 = \frac{3}{4}, \quad a_{2,1} = \frac{2}{3} \]公式：
\[y_{m+1} = y_m + \frac{h}{4} \left[ f(x_m, y_m) + 3 f\left(x_m + \frac{2h}{3}, y_m + \frac{2h}{3} f_m\right) \right] \]Butcher表：
\[\begin{array}{c|cc} 0 & & \\ \frac{2}{3} & \frac{2}{3} & \\ \hline & \frac{1}{4} & \frac{3}{4} \\ \end{array} \] -
改进Euler公式（取\(c_2 = 1\)）
\[b_1 = \frac{1}{2}, \quad b_2 = \frac{1}{2}, \quad a_{2,1} = 1 \]公式：
\[y_{m+1} = y_m + \frac{h}{2} \left[ f(x_m, y_m) + f(x_m + h, y_m + h f_m) \right] \]Butcher表：
\[\begin{array}{c|cc} 0 & & \\ 1 & 1 & \\ \hline & \frac{1}{2} & \frac{1}{2} \\ \end{array} \]
这正是我们上一节讲的预测-校正Euler方法！这说明改进Euler方法其实是二级二阶Runge-Kutta方法的一个特例。
四、常用的高阶Runge-Kutta方法
1. 三级三阶Runge-Kutta方法
对于\(s=3\)的情况，可以构造出三阶精度的方法。阶数条件共有6个方程，8个未知数，因此有两个自由参数。
(1) 三级三阶Kutta公式
Butcher表：
(2) 三级三阶Heun公式
Butcher表：
2. 四级四阶Runge-Kutta方法
这是工程中最常用的Runge-Kutta方法，通常称为经典四阶Runge-Kutta方法（RK4）。
(1) 经典四阶Runge-Kutta公式
Butcher表：
特点：
- 每步需要计算4次函数值
- 局部截断误差是\(O(h^5)\)，整体截断误差是\(O(h^4)\)
- 精度高，程序实现简单，是求解非刚性问题的首选方法
(2) 四级四阶Kutta公式
(3) 四级四阶Gill公式
Gill公式的优点：具有最小的舍入误差，特别适合长时间积分和高精度计算。
3. 更高阶的Runge-Kutta方法
- 五级五阶：不存在！这是一个重要的理论结果：5级显式Runge-Kutta方法最多只能达到4阶精度。
- 六级五阶：如Nyström公式、Lawson公式
- 七级六阶：如Butcher公式
阶数障碍：对于显式Runge-Kutta方法，当\(s \geq 5\)时，s级方法最多只能达到\(s-1\)阶精度。这是由Butcher在1963年证明的。
五、Butcher阵列的一般性质
对于任意显式Runge-Kutta方法，系数通常满足以下条件：
解释：
- 第一个条件\(\sum b_j = 1\)：保证方法是一阶精度的（一致性条件）
- 第二个条件\(c_i = \sum_{j=1}^{i-1} a_{i,j}\)：称为行和条件，保证了节点的一致性
六、核心知识点归纳总结
| 方法类型 | 级数\(s\) | 阶数\(p\) | 每步函数计算次数 | 局部截断误差 | 绝对稳定区间（实\(\lambda<0\)） | 典型公式 | 适用场景 |
|---|---|---|---|---|---|---|---|
| 显式Euler | 1 | 1 | 1 | \(O(h^2)\) | \(0 < h < -\frac{2}{\lambda}\) | \(y_{m+1}=y_m+hf_m\) | 精度要求不高的简单问题 |
| 二级二阶RK | 2 | 2 | 2 | \(O(h^3)\) | \(0 < h < -\frac{2}{\lambda}\) | 中点公式、Heun公式、改进Euler | 一般精度要求的问题 |
| 三级三阶RK | 3 | 3 | 3 | \(O(h^4)\) | \(0 < h < -\frac{2.51}{\lambda}\) | Kutta公式、Heun公式 | 中等精度要求的问题 |
| 四级四阶RK | 4 | 4 | 4 | \(O(h^5)\) | \(0 < h < -\frac{2.78}{\lambda}\) | 经典RK4、Gill公式 | 高精度要求的非刚性问题 |
| 六级五阶RK | 6 | 5 | 6 | \(O(h^6)\) | \(0 < h < -\frac{3.21}{\lambda}\) | Nyström公式、Lawson公式 | 高精度要求的长时间积分 |
七、学习建议
- 重点掌握经典四阶Runge-Kutta方法：这是工程中最常用的方法，必须熟练掌握其公式和编程实现。
- 理解阶数条件的推导：通过推导二级二阶方法的阶数条件，理解Runge-Kutta方法的构造原理。
- 熟悉Butcher表的使用：Butcher表是表示Runge-Kutta方法的标准方式，能够根据Butcher表写出对应的计算公式。
- 了解阶数障碍：知道为什么不存在五级五阶显式Runge-Kutta方法，这是数值分析中的一个重要理论结果。
- 通过编程实践加深理解：实现经典四阶Runge-Kutta方法，并用它求解不同的初值问题，观察步长对误差的影响。
1.3.2 隐式Runge-Kutta方法 完整讲解与深度证明
作为拥有多年教学经验的资深高级研究员，我将对这部分内容进行系统讲解，完整证明解的存在唯一性定理，揭示隐式Runge-Kutta方法的本质优势，并详细介绍工程中求解刚性问题最常用的Gauss-Legendre方法。
一、隐式Runge-Kutta方法的定义与分类
1. 一般形式
定义1.3.2 s级Runge-Kutta方法的一般形式为：
2. 分类（根据矩阵\(A=(a_{i,j})\)的结构）
Runge-Kutta方法的分类完全由内部权重矩阵\(A\)的结构决定：
| 方法类型 | 矩阵\(A\)的结构 | 特点 |
|---|---|---|
| 显式Runge-Kutta方法 | 严格下三角矩阵（\(a_{i,j}=0\)当\(i \leq j\)） | 每个\(k_i\)只依赖于前面的\(k_1, \dots, k_{i-1}\)，可以显式计算 |
| 对角隐式Runge-Kutta方法(DIERK) | 下三角矩阵且至少有一个对角元非零（\(a_{i,j}=0\)当\(i < j\)） | 每个\(k_i\)只依赖于\(k_1, \dots, k_i\)，可以逐个求解 |
| 单对角隐式Runge-Kutta方法(SDIRK) | 下三角矩阵且所有对角元相等（\(a_{i,i}=\gamma\)） | 所有\(k_i\)的求解具有相同的结构，便于编程实现 |
| 全隐式Runge-Kutta方法 | 一般矩阵（无零元素限制） | 所有\(k_i\)相互耦合，需要同时求解一个\(s\)维非线性方程组 |
3. Butcher表的统一形式
所有Runge-Kutta方法都可以用Butcher表统一表示：
显式方法的Butcher表：对角线及以上元素全为0
对角隐式方法的Butcher表：对角线以上元素全为0
全隐式方法的Butcher表：所有元素都可能非零
二、定理1.3.1（解的存在唯一性与迭代收敛性）完整证明
定理原文
定理1.3.1 设\(f(x,y)\)是连续的，且满足Lipschitz条件，\(L\)为相应的常数。如果
则方程(1.3.11)存在唯一解，并可通过迭代得到。进一步地，如果\(f(x,y)\)是\(p\)次连续可微函数，则\(k_i\)是\(h\)的\(p\)次连续可微函数。
完整证明过程
证明：我们将证明分为三个部分：存在唯一性、迭代收敛性和可微性。
1. 问题转化为不动点问题
记\(s\)维向量\(\boldsymbol{K} = (k_1, k_2, \dots, k_s)^T\)，定义算子\(\boldsymbol{F}: \mathbb{R}^s \to \mathbb{R}^s\)为：
其中
则方程(1.3.11)中的后\(s\)个方程可以写为不动点方程：
2. 证明\(\boldsymbol{F}\)是压缩映射
我们在\(\mathbb{R}^s\)上定义最大值范数：
现在计算\(\boldsymbol{F}\)的Lipschitz常数。对于任意两个向量\(\boldsymbol{K}_1 = (k_1^1, \dots, k_s^1)^T\)和\(\boldsymbol{K}_2 = (k_1^2, \dots, k_s^2)^T\)，有：
利用\(f\)关于\(y\)的Lipschitz条件：
得：
因此：
令
则：
根据定理条件(1.3.12)，\(\alpha < 1\)，因此\(\boldsymbol{F}\)是压缩映射。
3. 不动点的存在唯一性与迭代收敛性
根据压缩映射原理（Banach不动点定理），压缩映射\(\boldsymbol{F}\)在\(\mathbb{R}^s\)上存在唯一的不动点\(\boldsymbol{K}^*\)，即方程\(\boldsymbol{K} = \boldsymbol{F}(\boldsymbol{K})\)存在唯一解。
并且，不动点迭代
收敛于唯一解\(\boldsymbol{K}^*\)，收敛速度为线性收敛，收敛因子为\(\alpha\)。
分量形式的迭代公式为：
4. 可微性证明
将方程\(\boldsymbol{K} = \boldsymbol{F}(\boldsymbol{K})\)改写为：
我们需要证明\(\boldsymbol{K}(h)\)在\(h=0\)的邻域内是\(p\)次连续可微的。
当\(h=0\)时，\(\boldsymbol{\Phi}(0, \boldsymbol{K}) = \boldsymbol{K} - \boldsymbol{F}(\boldsymbol{K}) = \boldsymbol{K} - (f(x_m, y_m), \dots, f(x_m, y_m))^T\)，因此方程\(\boldsymbol{\Phi}(0, \boldsymbol{K}) = 0\)的解为：
现在计算\(\boldsymbol{\Phi}\)关于\(\boldsymbol{K}\)的Jacobi矩阵在\(h=0\)处的值：
当\(h=0\)时，\(\frac{\partial F_i}{\partial k_j}(0, \boldsymbol{K}) = 0\)对所有\(i,j\)成立，因此：
这是一个非奇异矩阵。
根据隐函数定理，如果\(\boldsymbol{\Phi}(h, \boldsymbol{K})\)是\(p\)次连续可微的，且\(\frac{\partial \boldsymbol{\Phi}}{\partial \boldsymbol{K}}(0, \boldsymbol{K}^*)\)非奇异，则存在\(h=0\)的一个邻域，在这个邻域内方程\(\boldsymbol{\Phi}(h, \boldsymbol{K}) = 0\)有唯一解\(\boldsymbol{K}(h)\)，且\(\boldsymbol{K}(h)\)是\(p\)次连续可微的。
由于\(f(x,y)\)是\(p\)次连续可微的，因此\(\boldsymbol{\Phi}(h, \boldsymbol{K})\)也是\(p\)次连续可微的。这就证明了\(k_i\)是\(h\)的\(p\)次连续可微函数。
证毕。
定理的重要意义
- 这个定理保证了隐式Runge-Kutta方法在步长\(h\)满足条件(1.3.12)时是适定的
- 它提供了一种求解隐式Runge-Kutta方法的标准方法：不动点迭代
- 可微性是进行误差分析和稳定性分析的基础
三、Gauss-Legendre方法：最常用的全隐式Runge-Kutta方法
隐式Runge-Kutta方法的最大优势是可以构造出任意高阶精度的方法，其中最著名的是Gauss-Legendre方法，它基于Gauss求积公式构造。
1. 基本思想
回顾积分形式的初值问题：
如果我们用s点Gauss求积公式近似计算右边的积分：
其中\(c_i\)是Gauss求积节点，\(b_i\)是Gauss求积权重。
Gauss求积公式的代数精度是\(2s-1\)，这意味着对于次数不超过\(2s-1\)的多项式，求积公式是精确的。对应的Runge-Kutta方法可以达到2s阶精度，这是显式Runge-Kutta方法无法比拟的（显式s级方法最多只能达到s阶精度）。
2. 参数确定方法
Gauss-Legendre方法的参数可以通过以下步骤确定：
步骤1：确定节点\(c_i\)
节点\(c_1, c_2, \dots, c_s\)是s次Legendre多项式\(P_s(2c - 1) = 0\)在区间\((0,1)\)内的根。
Legendre多项式的递推公式为：
步骤2：确定内部权重\(a_{i,j}\)
对于每个\(i = 1, 2, \dots, s\)，解线性方程组：
步骤3：确定外部权重\(b_i\)
解线性方程组：
3. 精度阶数
Butcher在1964年证明了：s级Gauss-Legendre方法的精度阶数是2s。
这是隐式Runge-Kutta方法最突出的优势：
- s=1：2阶精度
- s=2：4阶精度
- s=3：6阶精度
- s=4：8阶精度
- ...
而显式Runge-Kutta方法：
- s=1：1阶精度
- s=2：2阶精度
- s=3：3阶精度
- s=4：4阶精度
- s=5：4阶精度（阶数障碍）
四、典型的Gauss-Legendre方法
1. 隐式中点公式（s=1, p=2）
这是最简单的隐式Runge-Kutta方法，也是最常用的二阶隐式方法。
Butcher表：
计算公式：
特点：
- 2阶精度
- A-稳定（对任意步长\(h>0\)都绝对稳定）
- 每步需要解一个非线性方程
- 特别适合求解刚性问题
2. Hammer和Hollingsworth公式（s=2, p=4）
这是最常用的四阶隐式Runge-Kutta方法。
Butcher表：
计算公式：
特点：
- 4阶精度
- A-稳定
- 每步需要解一个二维非线性方程组
- 精度与经典四阶显式RK方法相同，但稳定性好得多
3. Kuntzmann和Butcher公式（s=3, p=6）
这是六阶隐式Runge-Kutta方法。
Butcher表：
特点：
- 6阶精度
- A-稳定
- 每步需要解一个三维非线性方程组
- 适合高精度要求的刚性问题
五、隐式Runge-Kutta方法的优缺点与适用场景
1. 核心优势
- 精度高：s级方法可以达到2s阶精度，远高于显式方法
- 稳定性好：所有Gauss-Legendre方法都是A-稳定的，对任意步长\(h>0\)都绝对稳定
- 适合刚性问题：这是隐式Runge-Kutta方法最主要的应用场景。对于刚性问题，显式方法必须使用非常小的步长才能保证稳定，而隐式方法可以使用大得多的步长，从而大大减少计算量
2. 主要缺点
- 计算量大：每步需要解一个s维非线性方程组，计算量远大于显式方法
- 编程复杂：需要实现非线性方程组的求解器（如牛顿法）
- 存储需求大：需要存储Jacobi矩阵和LU分解的结果
3. 显式与隐式Runge-Kutta方法对比
| 对比项 | 显式Runge-Kutta方法 | 隐式Runge-Kutta方法 |
|---|---|---|
| 精度 | s级最多s阶 | s级可达2s阶 |
| 稳定性 | 条件稳定，绝对稳定区域有限 | A-稳定，对任意步长都稳定 |
| 计算量 | 小，每步计算s次函数值 | 大，每步需要解s维非线性方程组 |
| 编程难度 | 简单 | 复杂，需要非线性求解器 |
| 适用场景 | 非刚性问题，精度要求一般 | 刚性问题，高精度要求 |
六、核心知识点归纳总结
| 知识点 | 核心内容 | 关键公式/结论 |
|---|---|---|
| 隐式RK定义 | 一般形式与分类 | 按矩阵\(A\)的结构分为显式、对角隐式、单对角隐式、全隐式 |
| 定理1.3.1 | 解的存在唯一性与迭代收敛性 | 收敛条件：\(hL \max_i \sum_j |a_{i,j}| < 1\) 证明方法：压缩映射原理 |
| Gauss-Legendre方法 | 基于Gauss求积构造的全隐式RK方法 | 节点：\(P_s(2c-1)=0\)的根 精度：2s阶 稳定性：A-稳定 |
| 典型方法 | 隐式中点公式(s=1,p=2) | \(y_{m+1}=y_m+hf(x_m+h/2,y_m+hk_1/2)\) |
| Hammer-Hollingsworth公式(s=2,p=4) | 四阶A-稳定方法 | |
| Kuntzmann-Butcher公式(s=3,p=6) | 六阶A-稳定方法 | |
| 优缺点 | 优势：高精度、A-稳定、适合刚性问题 缺点：计算量大、需要解非线性方程组 |
学习建议
- 重点掌握隐式中点公式：这是最简单也是最常用的隐式RK方法，必须熟练掌握其公式和求解方法。
- 理解压缩映射原理的应用：这是证明隐式方法解的存在唯一性的标准方法，也是数值分析中最基本的工具之一。
- 认识隐式方法的核心优势：A-稳定性是隐式方法的灵魂，也是它们能够求解刚性问题的根本原因。
- 通过编程实践加深理解：实现隐式中点公式，并用它求解一个刚性问题（如\(y' = -1000y\)），观察步长对稳定性的影响。
1.3.3 半隐式Runge-Kutta方法（Rosenbrock方法）完整讲解
作为拥有多年教学经验的资深高级研究员，我将对这部分内容进行系统讲解，揭示半隐式方法的本质思想，澄清"半隐式"的真正含义，并说明它在刚性问题求解中的核心优势。
一、半隐式方法的提出背景与核心思想
1. 现有方法的困境
在介绍半隐式方法之前，我们先回顾一下显式和隐式Runge-Kutta方法的优缺点：
| 方法类型 | 优点 | 缺点 |
|---|---|---|
| 显式RK | 计算简单，不需要解方程组，每步只需要计算s次函数值 | 稳定性差，绝对稳定区域有限，不适合求解刚性问题 |
| 隐式RK | 稳定性极好（A-稳定），精度高，适合求解刚性问题 | 每步需要解一个s维非线性方程组，计算量大，编程复杂 |
核心矛盾：对于刚性问题，显式方法步长太小导致计算量爆炸，隐式方法虽然可以用大步长，但每步的非线性迭代开销太大。
2. Rosenbrock的天才想法
1963年，德国数学家H. H. Rosenbrock提出了一个革命性的折中方案：在每个Runge-Kutta阶段进行一次牛顿迭代的线性化。
核心思想：
- 不直接求解非线性方程
- 而是对非线性函数\(f(x,y)\)在当前点进行泰勒展开，保留到一阶项（线性化）
- 这样就将每个阶段的非线性方程转化为线性方程
- 线性方程只需要求解一次，不需要迭代
本质：半隐式Runge-Kutta方法是显式计算流程 + 隐式线性阶段的结合体。它既保持了隐式方法的良好稳定性，又避免了非线性迭代的巨大开销。
二、半隐式Runge-Kutta方法的一般形式
1. 定义与公式
s级半隐式Runge-Kutta方法（Rosenbrock方法）的一般形式为：
其中\(k_i\)满足：
式中\(i = 2, 3, \dots, s\)，\(A(x,y) = f_y'(x,y)\)是\(f\)关于\(y\)的雅可比矩阵。
2. 关键参数解读
- \(a_{i,j}, \bar{a}_{i,j}\)：内部权重系数
- \(b_i\)：外部权重系数
- \(c_i, \bar{c}_i\)：节点参数
- \(\omega_i\)：线性化参数
- \(A(x,y) = f_y'(x,y)\)：雅可比矩阵，对于单个方程就是导数\(f_y(x,y)\)
3. "半隐式"的本质：显式计算流程
将(1.3.17)式整理一下，把含\(k_i\)的项移到左边：
因此：
关键观察：
- 计算\(k_i\)时，右边所有项都只依赖于前面已经计算出来的\(k_1, k_2, \dots, k_{i-1}\)
- 因此，我们可以按顺序\(k_1 \to k_2 \to \dots \to k_s\)逐个计算
- 每个\(k_i\)只需要求解一个线性方程组，不需要迭代
这就是"半隐式"的真正含义：每个阶段的方程是隐式的（线性），但整体计算流程是显式的。
4. 与显式RK的关系
当所有线性化参数\(\omega_i = 0\)时，(1.3.17)式退化为：
这正是显式Runge-Kutta方法的一般形式。因此，显式RK是半隐式RK当\(\omega_i = 0\)时的特例。
三、二级三阶Rosenbrock方法详解
1. 具体公式
对于\(s=2\)的情况，半隐式RK方法的公式为：
2. 单个方程的简化
对于单个常微分方程（不是方程组），雅可比矩阵\(A(x,y)\)退化为一个标量\(f_y(x,y)\)，因此矩阵求逆退化为标量求倒数：
重大优势：对于单个方程，半隐式RK方法完全不需要解任何方程组，只需要进行简单的算术运算！这是教材中特别强调的一点："对单个方程，它既保持隐式Runge-Kutta方法的好的稳定性，又不需要解非线性方程(组)"。
3. 二级三阶方法的系数
通过泰勒展开比较系数，可以构造出三阶精度的二级半隐式RK方法。其系数为：
4. 计算流程（单个方程）
以单个方程为例，二级三阶Rosenbrock方法的完整计算流程为：
- 计算当前点的函数值和导数值：\(f_m = f(x_m, y_m)\)，\(f_{y,m} = f_y(x_m, y_m)\)
- 计算\(k_1\)：\(k_1 = \frac{f_m}{1 - \omega_1 f_{y,m}}\)
- 计算中间点的函数值和导数值：\[x^* = x_m + c_2 h, \quad y^* = y_m + h a_{2,1} k_1 \]\[f^* = f(x^*, y^*), \quad f_y^* = f_y(x^*, y^*) \]
- 计算\(k_2\)：\(k_2 = \frac{f^*}{1 - \omega_2 f_y^*}\)
- 计算下一个点的解：\(y_{m+1} = y_m + h(b_1 k_1 + b_2 k_2)\)
可以看到，整个过程没有任何迭代，完全是显式计算！
四、半隐式Runge-Kutta方法的优缺点与适用场景
1. 核心优势
- 稳定性好：大多数Rosenbrock方法是L-稳定的（比A-稳定更强的稳定性），非常适合求解刚性问题
- 计算量小：不需要解非线性方程组，只需要解线性方程组（单个方程甚至不需要解方程组）
- 编程相对简单：比隐式RK方法容易实现得多
- 精度适中：s级方法可以达到s阶或s-1阶精度，对于大多数工程问题已经足够
2. 主要缺点
- 需要计算雅可比矩阵：这是半隐式方法最大的缺点。对于复杂的\(f(x,y)\)，计算雅可比矩阵可能比较麻烦
- 系数复杂：高阶Rosenbrock方法的系数非常复杂，而且没有统一的表达式
- 对于方程组需要解线性方程组：对于n维方程组，每个阶段需要解一个n阶线性方程组，虽然比非线性方程组简单，但还是比显式RK计算量大
3. 适用场景
半隐式Runge-Kutta方法是求解中等规模刚性问题的首选方法，特别适合以下情况：
- 单个或小规模刚性常微分方程组
- 对计算速度要求较高，同时需要良好稳定性的问题
- 雅可比矩阵容易计算的问题
在实际工程中，半隐式RK方法广泛应用于：
- 化学动力学模拟（反应速率常数相差悬殊）
- 电路模拟（不同时间尺度的电路元件）
- 控制系统仿真（刚性控制对象）
- 计算流体力学（对流扩散问题）
五、三种Runge-Kutta方法的全面对比
| 对比项 | 显式Runge-Kutta | 半隐式Runge-Kutta（Rosenbrock） | 隐式Runge-Kutta（Gauss-Legendre） |
|---|---|---|---|
| 计算流程 | 完全显式，按顺序计算\(k_i\) | 显式流程，每个阶段解线性方程组 | 隐式，需要解非线性方程组 |
| 迭代需求 | 不需要迭代 | 不需要迭代 | 需要牛顿迭代 |
| 雅可比矩阵 | 不需要 | 需要 | 需要（牛顿迭代） |
| 每步计算量 | 小，s次函数求值 | 中等，s次函数求值 + s次线性方程组求解 | 大，多次函数求值 + 多次线性方程组求解 |
| 稳定性 | 条件稳定，绝对稳定区域有限 | L-稳定，适合刚性问题 | A-稳定，最稳定 |
| 精度 | s级最多s阶 | s级可达s阶 | s级可达2s阶 |
| 编程难度 | 最简单 | 中等 | 最复杂 |
| 适用场景 | 非刚性问题 | 中等规模刚性问题 | 高精度刚性问题 |
六、扩展：W-方法（近似雅可比Rosenbrock方法）
为了克服需要计算精确雅可比矩阵的缺点，人们发展了W-方法（也称为"线化隐式RK方法"）。
基本思想：用一个近似雅可比矩阵\(W\)代替精确雅可比矩阵\(A(x,y)\)，并且在整个步长内保持不变。
优势：
- 不需要在每个阶段都重新计算雅可比矩阵
- 甚至可以用数值微分来近似雅可比矩阵
- 大大减少了计算量
代价：精度会略有降低，但稳定性通常保持不变。
W-方法是实际数值计算软件中最常用的刚性求解器之一，例如MATLAB中的ode23s和ode23t就是基于W-方法实现的。
七、核心知识点归纳总结
| 知识点 | 核心内容 | 关键公式/结论 |
|---|---|---|
| 提出背景 | 解决显式RK稳定性差、隐式RK计算量大的矛盾 | 刚性问题求解的折中方案 |
| 核心思想 | 每个阶段进行一次牛顿迭代的线性化，将非线性方程转化为线性方程 | 显式计算流程 + 隐式线性阶段 |
| 一般形式 | \(y_{m+1}=y_m+h\sum b_i k_i\) \(k_i=[I-\omega_i A]^{-1}f(\dots)\) |
|
| 半隐式本质 | 每个\(k_i\)只依赖于前面的\(k_1,\dots,k_{i-1}\)，可以按顺序计算 | 不需要迭代 |
| 单个方程简化 | 矩阵求逆退化为标量求倒数，完全不需要解方程组 | \(k_i = \frac{f(\dots)}{1 - \omega_i f_y(\dots)}\) |
| 典型方法 | 二级三阶Rosenbrock方法 | 系数见教材 |
| 优势 | L-稳定，适合刚性问题，计算量比隐式RK小 | |
| 缺点 | 需要计算雅可比矩阵 | |
| 适用场景 | 中等规模刚性问题 | 化学动力学、电路模拟、控制系统 |
学习建议
- 重点理解"半隐式"的本质：半隐式不是显式也不是隐式，而是显式流程与隐式线性阶段的结合。抓住"不需要迭代"这个核心特点。
- 掌握单个方程的计算流程：对于单个方程，半隐式RK方法非常简单，完全不需要解方程组，这是它的一大优势。
- 认识雅可比矩阵的作用：雅可比矩阵是半隐式方法的关键，它提供了关于解的局部变化率的信息，从而保证了方法的稳定性。
- 了解实际应用：半隐式RK方法是求解刚性问题的主流方法之一，很多数值计算软件都提供了相应的实现。
1.3.4 单步法的稳定性和收敛性 完整讲解与理论框架
作为拥有多年教学经验的资深高级研究员，我将把前面分散介绍的各种单步法（Euler、Runge-Kutta、半隐式RK）统一到一个通用理论框架下，系统讲解单步法的三个核心性质——相容性、稳定性、收敛性，并证明数值分析中最重要的定理之一：Lax等价定理。
一、单步法的一般形式与增量函数
1. 统一表示
所有单步法（无论显式还是隐式，无论一阶还是高阶）都可以写成如下统一形式：
其中\(\varphi(x_m, y_m, h)\)称为增量函数，它表示单位步长上的增量。
2. 常见方法的增量函数
我们可以将前面学过的所有单步法都表示为这种形式：
| 方法名称 | 增量函数\(\varphi(x,y,h)\) |
|---|---|
| 显式Euler方法 | \(\varphi(x,y,h) = f(x,y)\) |
| 隐式Euler方法 | \(\varphi(x,y,h) = f(x+h, y + h \varphi(x,y,h))\) |
| 改进Euler方法 | \(\varphi(x,y,h) = \frac{1}{2}\left[ f(x,y) + f(x+h, y + h f(x,y)) \right]\) |
| 二级二阶中点公式 | \(\varphi(x,y,h) = f\left(x + \frac{h}{2}, y + \frac{h}{2} f(x,y)\right)\) |
| 经典四阶RK方法 | \(\varphi(x,y,h) = \frac{1}{6}(k_1 + 2k_2 + 2k_3 + k_4)\) 其中\(k_1=f(x,y), k_2=f(x+h/2,y+hk_1/2), \dots\) |
意义：增量函数的引入将各种形式各异的单步法统一到了一个框架下，使得我们可以用一套通用的理论来分析它们的性质。
二、局部截断误差与精度阶数
1. 局部截断误差的一般定义
定义：假设前一步的计算值是精确的，即\(y_m = y(x_m)\)，则单步法的局部截断误差为：
2. 精度阶数的定义
定义1.3.3 如果\(p\)是使下式成立的最大整数：
则称单步法(1.3.19)是p阶精度的。
解读：
- 局部截断误差是\(O(h^{p+1})\)，意味着方法在单步计算中产生的误差与步长的\(p+1\)次方成正比
- 全局截断误差通常是\(O(h^p)\)，比局部误差低一阶
- \(p\)越大，方法的精度越高
三、相容性：方法与微分方程的一致性
1. 相容性的直观意义
一个数值方法要想逼近微分方程，当步长\(h \to 0\)时，它必须"收敛到"微分方程本身。
对于微分方程\(y' = f(x,y)\)，我们有：
而单步法可以写为：
因此，当\(h \to 0\)时，为了使数值方法逼近微分方程，必须有：
2. 相容性的严格定义
定义1.3.4 如果单步法的增量函数\(\varphi(x,y,h)\)满足：
则称单步法(1.3.19)与微分方程(1.1.1)相容。方程(1.3.20)称为相容性条件。
前提：\(\varphi(x,y,h)\)关于\(h\)在\(h=0\)处连续。
3. 相容性与精度的关系
定理：相容的单步法至少是一阶精度的。
证明：
将\(y(x+h)\)在\(x\)处展开为Taylor级数：
局部截断误差为：
如果方法相容，即\(\varphi(x,y,0) = f(x,y)\)，则：
因此：
这说明方法至少是一阶精度的。
结论：相容性是方法具有一阶精度的充要条件。
四、稳定性：初始误差的传播
1. 稳定性的一般定义
定义：称单步法(1.3.19)是稳定的，如果存在常数\(C\)和\(h_0\)，使得对于任意两个初始值\(y_0\)和\(z_0\)，以及任意\(0 < h < h_0\)，方法的解\(y_m\)和\(z_m\)满足：
2. 一般单步法的稳定性定理
定理1.3.2 如果对任意\((x,y) \in \Omega\)及\(0 < h \leq h_0\)，增量函数\(\varphi(x,y,h)\)关于\(y\)满足Lipschitz条件，即存在常数\(L > 0\)，使得：
对所有\((x,y_1), (x,y_2) \in \Omega\)成立，则单步法(1.3.19)是稳定的。
证明：与Euler方法的稳定性证明完全相同，只需将\(f(x,y)\)换成\(\varphi(x,y,h)\)即可。
考虑两个解：
两式相减，记误差\(e_m = y_m - z_m\)，得：
两边取绝对值，并利用Lipschitz条件：
递推得：
令\(C = e^{L(b-a)}\)，则：
证毕。
五、收敛性：数值解趋向于精确解
1. 收敛性的定义
定义1.3.5 称单步法(1.3.19)是收敛的，如果对于任意初值问题(1.1.1)-(1.1.2)和任意\(x \in [a,b]\)，有：
其中\(y(x)\)是初值问题的精确解。
解读：收敛性意味着当步长无限减小时，数值解在任意固定点\(x\)处趋向于精确解。
2. Lax等价定理（数值分析的基石）
定理1.3.3 在定理1.3.2的条件下，如果\(\varphi(x,y,h)\)关于\(x\)和\(h\)也满足Lipschitz条件，则单步法(1.3.19)的收敛性与相容性等价。
定理解读：
- 对于一个相容的单步法，收敛性等价于稳定性
- 换句话说：相容 + 稳定 ⇔ 收敛
这是数值分析中最重要的定理之一，它告诉我们：要判断一个单步法是否收敛，只需要验证两个条件：
- 方法是否与微分方程相容（\(\varphi(x,y,0)=f(x,y)\)）
- 方法是否稳定（增量函数满足Lipschitz条件）
只要这两个条件满足，方法就一定收敛。
3. 定理证明思路
证明分为两个方向：
- 充分性：相容 + 稳定 ⇒ 收敛
- 必要性：收敛 ⇒ 相容
充分性证明概要：
- 构造辅助微分方程：\(z' = g(x,z)\)，其中\(g(x,z) = \varphi(x,z,0)\)
- 证明数值解\(y_m\)收敛于辅助方程的解\(z(x)\)
- 如果方法相容，则\(g(x,z) = f(x,z)\)，因此辅助方程与原方程相同
- 因此数值解收敛于原方程的解\(y(x)\)
必要性证明概要：
- 假设数值解收敛于原方程的解\(y(x)\)
- 则数值解也收敛于辅助方程的解\(z(x)\)
- 因此\(y(x) = z(x)\)对所有\(x\)成立
- 求导得\(y'(x) = z'(x)\)，即\(f(x,y(x)) = g(x,y(x)) = \varphi(x,y(x),0)\)
- 由于初值\(y_0\)是任意的，因此对所有\((x,y)\)有\(\varphi(x,y,0) = f(x,y)\)，即方法相容
六、全局截断误差估计
1. 一般单步法的全局误差估计
定理1.3.4 在定理1.3.3的条件下，如果局部截断误差\(R_m\)满足：
则单步法的整体截断误差\(\varepsilon_m = y(x_m) - y_m\)满足估计式：
解读：
- 全局误差由两部分组成：
- 第一部分：\(e^{L(b-a)} |\varepsilon_0|\)，是初始误差的传播项
- 第二部分：\(h^p \frac{C}{L} \left( e^{L(b-a)} - 1 \right)\)，是局部截断误差的累积项
- 当初始误差\(\varepsilon_0 = 0\)时，全局误差是\(O(h^p)\)，与方法的阶数一致
2. 对Runge-Kutta方法的应用
前面介绍的所有Runge-Kutta方法都满足定理的条件：
- 当\(f(x,y)\)满足Lipschitz条件时，增量函数\(\varphi(x,y,h)\)也满足Lipschitz条件
- 因此，只要方法是相容的（所有RK方法都是相容的），它们就是收敛的
- p阶RK方法的全局截断误差是\(O(h^p)\)
七、核心理论框架总结
单步法的三个核心性质之间的逻辑关系可以用下图表示：
相容性(φ(x,y,0)=f(x,y)) ⇨ 至少一阶精度
↓
稳定性(φ关于y满足Lipschitz) ⇨ 初始误差不被无限放大
↓
Lax等价定理：相容 + 稳定 ⇔ 收敛
↓
p阶精度 + 收敛 ⇨ 全局误差O(h^p)
八、常见方法的三大性质验证
我们用这个通用理论来验证前面学过的各种方法：
| 方法名称 | 相容性 | 稳定性 | 收敛性 | 全局误差阶数 |
|---|---|---|---|---|
| 显式Euler | 是（\(\varphi=f\)） | 是（\(f\)满足Lipschitz） | 是 | \(O(h)\) |
| 隐式Euler | 是（\(\varphi=f(x+h,y+h\varphi)\)，\(\varphi(x,y,0)=f(x,y)\)） | 是 | 是 | \(O(h)\) |
| 改进Euler | 是 | 是 | 是 | \(O(h^2)\) |
| 经典RK4 | 是 | 是 | 是 | \(O(h^4)\) |
| 二级三阶Rosenbrock | 是 | 是 | 是 | \(O(h^3)\) |
九、学习建议
- 深刻理解Lax等价定理：这是数值分析的基石，它将三个看似独立的概念联系在了一起，是判断数值方法好坏的根本标准。
- 掌握增量函数的思想：增量函数是统一所有单步法的关键，能够将具体方法表示为增量函数形式是理解通用理论的前提。
- 区分局部误差和全局误差：局部误差是单步误差，全局误差是累积误差，全局误差的阶数比局部误差低一阶。
- 理解收敛性的本质：收敛性是当\(h \to 0\)时的极限性质，它保证了数值方法在理论上的正确性。
1.4 线性多步法 入门讲解与核心概念
作为拥有多年教学经验的资深高级研究员，我将系统讲解线性多步法的基本思想、一般形式和核心概念，对比它与单步法的本质区别，并揭示其在大规模数值计算中的独特优势。
一、多步法的基本思想：利用历史信息提高效率
1. 单步法的局限性
前面介绍的所有单步法（Euler、Runge-Kutta、半隐式RK）都有一个共同特点：计算\(y_{m+1}\)时只用到前一步的信息\(y_m\)。
为了提高精度，单步法需要在每一步内计算多个点的函数值（例如经典RK4每步计算4次\(f\)）。这对于精度要求高的问题，计算量会急剧增加。
2. 多步法的天才想法
多步法的核心思想非常朴素：既然我们已经计算出了前面\(k\)步的结果\(y_m, y_{m-1}, \dots, y_{m-k+1}\)，为什么不把这些历史信息利用起来呢？
通过合理组合前面\(k\)步的解\(y_i\)和导数值\(f_i = f(x_i, y_i)\)，我们可以构造出精度更高、每步只需要计算一次函数值的数值方法。这就是多步法的本质。
3. 多步法的核心优势
- 计算效率高：对于p阶精度的方法，线性多步法每步只需要计算1次\(f\)，而p阶Runge-Kutta方法每步需要计算p次\(f\)
- 存储需求小：只需要存储前面\(k\)步的\(y\)和\(f\)值
- 易于编程实现：公式形式简单，循环结构清晰
二、多步法的一般形式
1. 一般定义
k步多步法的一般形式为：
2. 参数含义
- \(k\)：步数，表示方法用到了前面\(k\)个点的信息
- \(\alpha_i\)：解的系数，是给定的实数
- \(\Phi\)：增量函数，类似于单步法的增量函数，但它依赖于前面\(k+1\)个点的\(x\)和\(f\)值
- \(f_i = f(x_i, y_i)\)：第\(i\)步的导数值
3. 显式与隐式多步法
- 显式多步法：如果增量函数\(\Phi\)不依赖于未知的\(f_{m+1}\)（即不依赖于\(y_{m+1}\)），则方法是显式的
- 隐式多步法：如果增量函数\(\Phi\)依赖于\(f_{m+1}\)，则方法是隐式的，每步需要解一个关于\(y_{m+1}\)的方程
三、局部截断误差与精度阶数
1. 局部截断误差的定义
与单步法类似，多步法的局部截断误差定义为：假设前面所有步的计算值都是精确的，即\(y_i = y(x_i)\)对所有\(i \leq m\)成立，则用多步法计算\(y_{m+1}\)产生的误差：
2. 精度阶数的定义
如果\(p\)是使下式成立的最大整数：
则称该多步法是p阶精度的。
解读：
- 局部截断误差是\(O(h^{p+1})\)，意味着单步误差与步长的\(p+1\)次方成正比
- 全局截断误差通常是\(O(h^p)\)，与单步法一致
四、线性多步法：最常用的多步法类型
1. 定义
如果多步法的增量函数\(\Phi\)是\(f\)的线性组合，即：
则该多步法称为线性多步法。
k步线性多步法的一般形式为：
其中\(\alpha_i\)和\(\beta_i\)是常数，且满足\(\alpha_k^2 + \beta_k^2 \neq 0\)（保证方法确实是k步的）。
2. 显式与隐式线性多步法
- 显式线性多步法：\(\beta_0 = 0\)，此时公式右边不包含未知的\(f_{m+1}\)，可以直接计算\(y_{m+1}\)
- 隐式线性多步法：\(\beta_0 \neq 0\)，此时公式右边包含未知的\(f_{m+1}\)，需要解关于\(y_{m+1}\)的方程
3. 线性多步法的特点
- 线性结构：\(y_{m+1}\)是前面\(k\)个\(y\)和\(k+1\)个\(f\)的线性组合
- 非自启动：k步线性多步法需要前面\(k\)个初始值\(y_0, y_1, \dots, y_{k-1}\)才能开始计算，而通常只有\(y_0\)是已知的，因此需要用单步法（如RK4）来计算初始的\(k-1\)个值
- 每步计算量小：显式线性多步法每步只需要计算1次\(f\)，隐式线性多步法每步通常只需要1-2次迭代
五、常见线性多步法举例
为了让大家有个直观的认识，这里给出几个最常用的线性多步法的例子：
1. 显式线性多步法：Adams-Bashforth方法
- 2步2阶Adams-Bashforth方法：\[y_{m+1} = y_m + \frac{h}{2}(3f_m - f_{m-1}) \]
- 3步3阶Adams-Bashforth方法：\[y_{m+1} = y_m + \frac{h}{12}(23f_m - 16f_{m-1} + 5f_{m-2}) \]
- 4步4阶Adams-Bashforth方法：\[y_{m+1} = y_m + \frac{h}{24}(55f_m - 59f_{m-1} + 37f_{m-2} - 9f_{m-3}) \]
2. 隐式线性多步法：Adams-Moulton方法
- 1步2阶Adams-Moulton方法（梯形方法）：\[y_{m+1} = y_m + \frac{h}{2}(f_{m+1} + f_m) \]
- 2步3阶Adams-Moulton方法：\[y_{m+1} = y_m + \frac{h}{12}(5f_{m+1} + 8f_m - f_{m-1}) \]
- 3步4阶Adams-Moulton方法：\[y_{m+1} = y_m + \frac{h}{24}(9f_{m+1} + 19f_m - 5f_{m-1} + f_{m-2}) \]
六、单步法与线性多步法的全面对比
| 对比项 | 单步法（如RK4） | 线性多步法（如4阶Adams） |
|---|---|---|
| 用到的信息 | 只用到前一步\(y_m\) | 用到前面\(k\)步的\(y\)和\(f\) |
| 每步函数计算次数 | p阶方法需要p次 | p阶方法只需要1次 |
| 自启动性 | 是，只需要初始值\(y_0\) | 否，需要k个初始值 |
| 变步长难度 | 容易，每步可以独立调整步长 | 较难，变步长需要重新计算历史值 |
| 编程复杂度 | 中等 | 简单 |
| 计算效率（高精度） | 低 | 高 |
| 适用场景 | 小规模问题、变步长需求高的问题 | 大规模问题、精度要求高的问题 |
七、核心知识点归纳总结
| 知识点 | 核心内容 | 关键公式/结论 |
|---|---|---|
| 基本思想 | 利用前面k步的历史信息提高计算效率 | |
| 多步法一般形式 | \(y_{m+1} = \sum_{i=1}^k \alpha_i y_{m-i+1} + h \Phi(\dots)\) | |
| 显式vs隐式 | 显式：\(\Phi\)不依赖\(y_{m+1}\) 隐式：\(\Phi\)依赖\(y_{m+1}\) |
|
| 局部截断误差 | 假设前面所有步精确时的单步误差 | \(R = O(h^{p+1})\) |
| 精度阶数 | 局部截断误差的阶数减1 | p阶方法全局误差\(O(h^p)\) |
| 线性多步法 | \(\Phi\)是\(f\)的线性组合 | \(y_{m+1} = \sum_{i=1}^k \alpha_i y_{m-i+1} + h \sum_{i=0}^k \beta_i f_{m-i+1}\) |
| 核心优势 | 每步只需要计算1次\(f\)，计算效率高 | |
| 主要缺点 | 非自启动，变步长较难 |
学习建议
- 重点理解多步法的基本思想：利用历史信息是多步法与单步法最本质的区别，也是其高效率的来源。
- 掌握线性多步法的一般形式：能够识别显式和隐式线性多步法，理解每个系数的含义。
- 认识线性多步法的优缺点：知道什么时候应该用多步法，什么时候应该用单步法。
- 了解常见的线性多步法：记住4阶Adams-Bashforth和Adams-Moulton方法的公式，它们是工程中最常用的线性多步法。
1.4.1 Adams外插法（Adams-Bashforth方法）完整推导与讲解
作为拥有多年教学经验的资深高级研究员，我将系统讲解Adams外插法的构造原理、完整推导过程、误差分析和常用公式，揭示其作为最经典显式线性多步法的核心优势与特点。
一、Adams方法的基本思想
1. 出发点：积分形式的初值问题
所有Adams方法（包括外插法和内插法）都源于初值问题的积分形式：
这是一个精确等式，没有任何近似。
2. 核心思路：用插值多项式近似被积函数
直接计算积分\(\int_{x_m}^{x_{m+1}} f(x, y(x)) dx\)是困难的，因为我们不知道\(y(x)\)的精确表达式。Adams方法的核心思想是：
- 利用已经计算出的前面\(k\)个点的导数值\(f_m, f_{m-1}, \dots, f_{m-k+1}\)
- 构造一个\(k-1\)次插值多项式\(L_{m,k-1}(x)\)来近似\(f(x, y(x))\)
- 然后对插值多项式积分，得到\(y_{m+1}\)的近似值
3. 外插法的定义
如果我们选取的插值节点是\(x_m, x_{m-1}, \dots, x_{m-k+1}\)，这些点都在积分区间\([x_m, x_{m+1}]\)的左侧，那么我们是在外插积分区间内的函数值，因此这种方法称为Adams外插法（也叫Adams显式方法）。
特点：插值多项式不依赖于未知的\(f_{m+1}\)，因此得到的数值方法是显式的，可以直接计算\(y_{m+1}\)。
二、基于Newton向后插值的公式推导
1. Newton向后插值公式
由于插值节点是等距的（步长为\(h\)），且是从\(x_m\)往回数的，使用Newton向后插值公式最为方便。
令\(x = x_m + \tau h\)，其中\(\tau \in [0,1]\)（对应积分区间\([x_m, x_{m+1}]\)）。Newton向后插值公式为：
其中\(\nabla^j f_m\)表示\(j\)阶向后差分，定义为：
- 0阶向后差分：\(\nabla^0 f_m = f_m\)
- 1阶向后差分：\(\nabla f_m = f_m - f_{m-1}\)
- 2阶向后差分：\(\nabla^2 f_m = \nabla f_m - \nabla f_{m-1} = f_m - 2f_{m-1} + f_{m-2}\)
- 一般地：\(\nabla^j f_m = \sum_{i=0}^j (-1)^i \binom{j}{i} f_{m-i}\)
2. 二项式系数的推广
为了简化表达式，我们引入广义二项式系数：
利用广义二项式系数，Newton向后插值公式可以简洁地表示为：
3. 积分得到差分形式的Adams公式
将插值多项式代入积分式(1.4.4)，并舍去余项，得到：
做变量替换\(x = x_m + \tau h\)，则\(dx = h d\tau\)，积分上下限从\(x \in [x_m, x_{m+1}]\)变为\(\tau \in [0,1]\)，因此：
定义系数：
则得到差分形式的k步Adams外插公式：
三、系数\(a_j\)的计算与递推公式
直接通过积分计算\(a_j\)比较麻烦，我们可以通过母函数方法推导出\(a_j\)的递推公式。
1. 母函数的定义
定义\(a_j\)的母函数为：
根据\(a_j\)的定义：
整理得：
2. Taylor级数展开与递推公式
我们知道以下两个Taylor级数展开式：
将这两个展开式代入(1.4.12)式，得到：
比较等式两边\(\tau^j\)的系数，得到\(a_j\)的递推公式：
3. 前几个\(a_j\)的值
利用递推公式，我们可以很容易地计算出前几个\(a_j\)的值：
| \(j\) | 0 | 1 | 2 | 3 | 4 | 5 |
|---|---|---|---|---|---|---|
| \(a_j\) | 1 | \(\frac{1}{2}\) | \(\frac{5}{12}\) | \(\frac{3}{8}\) | \(\frac{251}{720}\) | \(\frac{95}{288}\) |
四、局部截断误差分析
1. 插值余项
Newton向后插值的余项为：
2. 局部截断误差的推导
Adams外插法的局部截断误差是插值余项的积分：
代入余项表达式，做变量替换\(x = x_m + \tau h\)，得：
3. 结论
k步Adams外插法的局部截断误差是\(O(h^{k+1})\)，因此它是k阶精度的方法。
这是一个非常重要的结论：k步Adams外插法具有k阶精度。例如：
- 1步Adams外插法（Euler方法）：1阶精度
- 2步Adams外插法：2阶精度
- 3步Adams外插法：3阶精度
- 4步Adams外插法：4阶精度
五、函数值形式的Adams外插公式
在实际计算中，使用函数值比使用向后差分更方便。我们可以利用向后差分与函数值的关系：
将其代入差分形式的Adams公式(1.4.10)，得到函数值形式的k步Adams外插公式：
其中系数\(b_{k,i}\)由下式确定：
六、常用的Adams外插公式
利用上面的公式，我们可以计算出常用的1到4步Adams外插公式：
1. 1步1阶Adams外插法（显式Euler方法）
局部截断误差：\(R_m = \frac{h^2}{2} y''(\xi)\)
2. 2步2阶Adams外插法
局部截断误差：\(R_m = \frac{5h^3}{12} y'''(\xi)\)
3. 3步3阶Adams外插法
局部截断误差：\(R_m = \frac{3h^4}{8} y^{(4)}(\xi)\)
4. 4步4阶Adams外插法（最常用）
局部截断误差：\(R_m = \frac{251h^5}{720} y^{(5)}(\xi)\)
七、Adams外插法的重要性质
1. 相容性
所有Adams外插法都满足相容性条件：
这可以通过将\(h \to 0\)时公式退化为微分方程\(y' = f(x,y)\)来验证。
2. 自启动性
k步Adams外插法不是自启动的，它需要前面\(k\)个初始值\(y_0, y_1, \dots, y_{k-1}\)才能开始计算。通常我们用单步法（如经典四阶Runge-Kutta方法）来计算初始的\(k-1\)个值。
3. 稳定性
Adams外插法是条件稳定的，其绝对稳定区间随着步数的增加而减小。例如：
- 1步（Euler）：\((-2, 0)\)
- 2步：\((-1, 0)\)
- 3步：\((-0.55, 0)\)
- 4步：\((-0.3, 0)\)
这是显式线性多步法的共同特点：步数越多，精度越高，但绝对稳定区间越小。
八、核心知识点归纳总结
| 方法名称 | 步数k | 阶数p | 计算公式 | 局部截断误差主项 |
|---|---|---|---|---|
| 显式Euler | 1 | 1 | \(y_{m+1}=y_m+hf_m\) | \(\frac{h^2}{2}y''(\xi)\) |
| 2步Adams-Bashforth | 2 | 2 | \(y_{m+1}=y_m+\frac{h}{2}(3f_m-f_{m-1})\) | \(\frac{5h^3}{12}y'''(\xi)\) |
| 3步Adams-Bashforth | 3 | 3 | \(y_{m+1}=y_m+\frac{h}{12}(23f_m-16f_{m-1}+5f_{m-2})\) | \(\frac{3h^4}{8}y^{(4)}(\xi)\) |
| 4步Adams-Bashforth | 4 | 4 | \(y_{m+1}=y_m+\frac{h}{24}(55f_m-59f_{m-1}+37f_{m-2}-9f_{m-3})\) | \(\frac{251h^5}{720}y^{(5)}(\xi)\) |
学习建议
- 重点掌握4步4阶Adams外插公式：这是工程中最常用的显式线性多步法，必须熟记其公式和局部截断误差。
- 理解插值积分的构造思想：Adams方法是数值积分与插值多项式结合的典范，理解这一思想有助于掌握其他线性多步法。
- 注意非自启动性：使用多步法时，必须先用单步法计算初始值，这是实际编程中容易忽略的一点。
- 了解稳定性特点：显式Adams方法的绝对稳定区间较小，不适合求解刚性问题。
1.4.2 Adams内插法（Adams-Moulton方法）完整讲解
作为拥有60多年教学经验的资深高级研究员，我将系统讲解Adams内插法的构造原理、完整推导过程、误差分析和常用公式，并对比它与Adams外插法的本质区别，最后介绍工程中最常用的Adams预测-校正方法。
一、Adams内插法的基本思想与核心优势
1. 插值节点的选择：内插vs外插
Adams外插法和内插法都源于初值问题的积分形式：
两者的唯一区别在于插值节点的选择：
- Adams外插法：选取\(x_m, x_{m-1}, \dots, x_{m-k+1}\)共\(k\)个点作为插值节点，这些点都在积分区间\([x_m, x_{m+1}]\)的左侧，因此是外插
- Adams内插法：选取\(x_{m+1}, x_m, \dots, x_{m-k+1}\)共\(k+1\)个点作为插值节点，这些点包含了积分区间的右端点\(x_{m+1}\)，因此是内插
2. 内插法的核心优势
根据插值理论，同样次数的内插公式精度远高于外插公式。具体来说：
- k步Adams外插法使用\(k\)个点，构造\(k-1\)次多项式，精度为\(k\)阶
- k步Adams内插法使用\(k+1\)个点，构造\(k\)次多项式，精度为\(k+1\)阶
此外，Adams内插法的数值稳定性也远好于外插法，绝对稳定区间更大，更适合长时间积分。
3. 隐式特性
由于插值节点包含了未知的\(x_{m+1}\)，Adams内插法的公式中会出现未知的\(f_{m+1} = f(x_{m+1}, y_{m+1})\)，因此它是隐式方法，每步需要解一个关于\(y_{m+1}\)的方程。
二、Adams内插法的完整推导
1. 插值多项式的构造
我们选取\(k+1\)个等距节点\(x_{m+1}, x_m, \dots, x_{m-k+1}\)，构造一个\(k\)次Newton向后插值多项式\(L_{m,k}^{(1)}(x)\)来近似\(f(x, y(x))\)。
由于插值节点是从\(x_{m+1}\)往回数的，我们以\(x_{m+1}\)为基点进行向后插值。令：
其中\(\tau = -1\)对应\(x = x_m\)，\(\tau = 0\)对应\(x = x_{m+1}\)。
Newton向后插值公式为：
其中\(\nabla^j f_{m+1}\)是在\(x_{m+1}\)处的\(j\)阶向后差分。
2. 积分得到差分形式的公式
将插值多项式代入积分式，舍去余项，得到：
做变量替换\(x = x_{m+1} + \tau h\)，则\(dx = h d\tau\)，积分上下限从\(x \in [x_m, x_{m+1}]\)变为\(\tau \in [-1, 0]\)，因此：
定义内插法的系数：
则得到差分形式的k步Adams内插公式：
3. 系数\(a_j^*\)的递推公式
与外插法类似，我们可以通过母函数方法推导出\(a_j^*\)的递推公式。
定义\(a_j^*\)的母函数为：
根据\(a_j^*\)的定义，可以推导出：
结合Taylor级数展开：
比较系数，得到\(a_j^*\)的递推公式：
4. 前几个\(a_j^*\)的值
利用递推公式，计算出前几个\(a_j^*\)的值：
| \(j\) | 0 | 1 | 2 | 3 | 4 | 5 |
|---|---|---|---|---|---|---|
| \(a_j^*\) | 1 | \(-\frac{1}{2}\) | \(-\frac{1}{12}\) | \(-\frac{1}{24}\) | \(-\frac{19}{720}\) | \(-\frac{3}{160}\) |
5. 函数值形式的Adams内插公式
利用向后差分与函数值的关系：
将其代入差分形式的公式，得到函数值形式的k步Adams内插公式：
其中系数\(b_{k,i}^*\)由下式确定：
三、局部截断误差分析
1. 插值余项
k次Newton向后插值的余项为：
2. 局部截断误差的推导
Adams内插法的局部截断误差是插值余项的积分：
代入余项表达式，做变量替换\(x = x_{m+1} + \tau h\)，得：
3. 关键结论
k步Adams内插法的局部截断误差是\(O(h^{k+2})\)，因此它是k+1阶精度的方法。
这是Adams内插法最突出的优势：同步数下，内插法比外插法精度高一阶。例如：
- 1步内插法：2阶精度（外插法1阶）
- 2步内插法：3阶精度（外插法2阶）
- 3步内插法：4阶精度（外插法3阶）
- 4步内插法：5阶精度（外插法4阶）
四、常用的Adams内插公式
1. 1步2阶Adams内插法（梯形公式）
局部截断误差：\(R_m = -\frac{h^3}{12} y'''(\xi)\)
说明：这正是我们之前学过的梯形方法！这说明梯形方法其实是1步2阶Adams内插法。
2. 2步3阶Adams内插法
局部截断误差：\(R_m = -\frac{h^4}{24} y^{(4)}(\xi)\)
3. 3步4阶Adams内插法（最常用）
局部截断误差：\(R_m = -\frac{19h^5}{720} y^{(5)}(\xi)\)
4. 4步5阶Adams内插法
局部截断误差：\(R_m = -\frac{3h^6}{160} y^{(6)}(\xi)\)
五、Adams内插法的重要性质
1. 隐式特性与迭代求解
Adams内插法是隐式方法，公式右边包含未知的\(f_{m+1} = f(x_{m+1}, y_{m+1})\)，因此需要通过迭代求解。
迭代公式：
其中\(f_{m+1}^{[\nu]} = f(x_{m+1}, y_{m+1}^{[\nu]})\)，而\(f_m, f_{m-1}, \dots\)都是已知的。
初始猜测：通常用Adams外插法的结果作为初始猜测\(y_{m+1}^{[0]}\)。由于外插法已经有k阶精度，通常只需要1-2次迭代就可以达到内插法的k+1阶精度。
2. 稳定性优势
Adams内插法的绝对稳定区间远大于同阶的外插法：
| 方法 | 步数 | 阶数 | 绝对稳定区间（实\(\lambda<0\)） |
|---|---|---|---|
| Adams外插 | 4 | 4 | \((-0.3, 0)\) |
| Adams内插 | 3 | 4 | \((-3.0, 0)\) |
可以看到，3步4阶内插法的绝对稳定区间是4步4阶外插法的10倍！这使得内插法可以使用更大的步长，从而大大提高计算效率。
3. 与外插法的系数关系
外插法系数\(a_j\)和内插法系数\(a_j^*\)之间存在一个非常简洁的关系：
或者等价地：
这个关系可以通过母函数\(G^*(\tau) = (1 - \tau) G(\tau)\)直接推导出来。
六、Adams预测-校正方法（ABM方法）
纯隐式的Adams内插法需要迭代求解，计算量较大。在实际工程中，最常用的是Adams预测-校正方法，它结合了外插法的显式特性和内插法的高精度与稳定性。
1. 基本思想
- 预测步：用k步Adams外插法计算一个预测值\(\bar{y}_{m+1}\)
- 计算预测导数值：\(\bar{f}_{m+1} = f(x_{m+1}, \bar{y}_{m+1})\)
- 校正步：用k步Adams内插法，将预测值\(\bar{f}_{m+1}\)代入，计算校正值\(y_{m+1}\)
2. 4阶Adams预测-校正公式（最常用）
这是工程中最常用的线性多步法，由4步4阶Adams外插法（AB4）和3步4阶Adams内插法（AM4）组成：
预测步（AB4）：
计算预测导数：
校正步（AM4）：
特点：
- 每步只需要计算1次新的函数值（预测步的\(\bar{f}_{m+1}\)）
- 整体精度为4阶，与经典RK4相同
- 计算效率远高于RK4（RK4每步需要计算4次函数值）
- 稳定性远好于显式RK方法
七、Adams外插法与内插法全面对比
| 对比项 | Adams外插法（AB） | Adams内插法（AM） | Adams预测-校正（ABM） |
|---|---|---|---|
| 插值节点 | \(x_m, \dots, x_{m-k+1}\)（k个点） | \(x_{m+1}, \dots, x_{m-k+1}\)（k+1个点） | 外插预测+内插校正 |
| 显/隐式 | 显式 | 隐式 | 显式（一次校正） |
| 步数k | k步k阶 | k步k+1阶 | 4步4阶（最常用） |
| 每步函数计算次数 | 1次 | 迭代多次 | 1次 |
| 局部截断误差 | \(O(h^{k+1})\) | \(O(h^{k+2})\) | \(O(h^5)\)（4阶） |
| 绝对稳定区间 | 小（如4步：(-0.3,0)） | 大（如3步：(-3.0,0)） | 中等（约(-0.7,0)） |
| 自启动性 | 否 | 否 | 否 |
| 编程复杂度 | 简单 | 中等 | 简单 |
| 计算效率 | 高 | 中 | 最高 |
| 适用场景 | 非刚性问题，精度要求一般 | 高精度刚性问题 | 大规模非刚性问题 |
八、核心知识点归纳总结
| 方法名称 | 步数 | 阶数 | 计算公式 | 局部截断误差主项 |
|---|---|---|---|---|
| 梯形公式（AM1） | 1 | 2 | \(y_{m+1}=y_m+\frac{h}{2}(f_{m+1}+f_m)\) | \(-\frac{h^3}{12}y'''(\xi)\) |
| 2步AM | 2 | 3 | \(y_{m+1}=y_m+\frac{h}{12}(5f_{m+1}+8f_m-f_{m-1})\) | \(-\frac{h^4}{24}y^{(4)}(\xi)\) |
| 3步AM | 3 | 4 | \(y_{m+1}=y_m+\frac{h}{24}(9f_{m+1}+19f_m-5f_{m-1}+f_{m-2})\) | \(-\frac{19h^5}{720}y^{(5)}(\xi)\) |
| 4阶ABM | 4 | 4 | AB4预测+AM4校正 | \(O(h^5)\) |
学习建议
- 重点掌握4阶Adams预测-校正方法：这是工程中最常用的线性多步法，兼顾了精度、稳定性和计算效率。
- 理解内插法精度更高的原因：插值节点包含积分区间内的点，多项式次数更高，误差更小。
- 认识隐式方法的稳定性优势：内插法的绝对稳定区间远大于外插法，这是它最核心的优势之一。
- 掌握预测-校正的思想：用显式方法预测，隐式方法校正，是解决隐式方法迭代问题的标准技巧。
1.4.3 一般线性多步公式 完整讲解与待定系数法
作为拥有60多年教学经验的资深高级研究员，我将系统讲解一般线性多步公式的构造方法——待定系数法，这是构造任意阶线性多步法的通用工具。我将完整推导阶数条件，并通过两步法的例子详细展示待定系数法的应用，最后介绍著名的Milne方法。
一、线性多步法的一般形式
1. 最一般形式
所有线性多步法都可以写成如下统一形式：
其中\(f_{m+j} = f(x_{m+j}, y_{m+j})\)，\(\alpha_j\)和\(\beta_j\)是常数，且满足：
- \(\alpha_k = 1\)：保证方法是k步的，且当\(h\)充分小时隐式方程有解
- \(|\alpha_0| + |\beta_0| \neq 0\)：保证方法确实用到了前面的信息
2. 显式与隐式的判断
- 显式线性多步法：如果\(\beta_k = 0\)，则公式右边不包含未知的\(f_{m+k}\)，可以直接计算\(y_{m+k}\)
- 隐式线性多步法：如果\(\beta_k \neq 0\)，则公式右边包含未知的\(f_{m+k}\)，需要解关于\(y_{m+k}\)的方程
3. 与之前方法的关系
这个一般形式包含了我们之前学过的所有方法：
- 当\(k=1, \alpha_1=1, \alpha_0=-1, \beta_1=0, \beta_0=1\)时，就是显式Euler方法
- 当\(k=1, \alpha_1=1, \alpha_0=-1, \beta_1=1, \beta_0=0\)时，就是隐式Euler方法
- 当\(k=1, \alpha_1=1, \alpha_0=-1, \beta_1=\beta_0=1/2\)时，就是梯形方法
- 所有Adams外插法和内插法都是这个一般形式的特例
二、待定系数法：构造线性多步法的通用方法
插值法虽然直观，但只能构造基于积分的Adams类方法。待定系数法是一种更通用、更灵活的构造方法，可以构造出任意形式的线性多步法。
1. 基本思想
待定系数法的核心思想是：通过Taylor展开比较系数，确定线性多步法的系数\(\alpha_j\)和\(\beta_j\)，使得方法达到尽可能高的精度。
2. 线性差分算子的定义
为了系统地推导阶数条件，我们定义线性差分算子：
这个算子的意义是：将精确解\(y(x)\)代入线性多步法公式，得到的残差。如果方法是p阶精度的，那么这个残差应该是\(O(h^{p+1})\)。
3. Taylor展开与阶数条件
将\(y(x + jh)\)和\(y'(x + jh)\)在点\(x\)处展开为Taylor级数：
将这两个展开式代入算子\(L[y(x); h]\)，整理得：
其中系数\(c_i\)由下式确定：
4. 精度阶数的定义
如果存在最大的整数\(p\)，使得：
则称该线性多步法是p阶精度的。此时，局部截断误差为：
其中\(c_{p+1}\)称为局部截断误差主项系数。
5. 相容性条件
当\(p \geq 1\)时，必有\(c_0 = 0\)和\(c_1 = 0\)，即：
这两个条件称为相容性条件，是线性多步法具有至少一阶精度的必要条件。
三、待定系数法的应用：两步法的构造
我们以两步法（\(k=2\)）为例，详细展示待定系数法的应用过程。
1. 设定参数
对于两步法，\(k=2\)，且\(\alpha_2 = 1\)。我们设\(\alpha_0 = \alpha\)（自由参数），则需要确定的系数为\(\alpha_1, \beta_0, \beta_1, \beta_2\)，共4个未知数。
2. 建立阶数条件
为了使方法达到尽可能高的精度，我们令\(c_0 = c_1 = c_2 = c_3 = 0\)，得到4个方程：
3. 解方程组
解这个线性方程组，得到：
4. 得到一般两步法公式
将这些系数代入一般形式(1.4.31)，得到含一个自由参数\(\alpha\)的两步法公式：
5. 精度分析
我们可以进一步计算\(c_4\)和\(c_5\)：
由此得到精度结论：
- 当\(\alpha \neq -1\)时，\(c_4 \neq 0\)，因此方法是三阶精度的
- 当\(\alpha = -1\)时，\(c_4 = 0\)，\(c_5 \neq 0\)，因此方法是四阶精度的
四、Milne方法：著名的四阶两步法
当\(\alpha = -1\)时，我们得到一个非常著名的四阶两步法——Milne方法。
1. Milne方法的公式
将\(\alpha = -1\)代入(1.4.36)，得到：
2. 局部截断误差
Milne方法的局部截断误差主项为：
3. 方法特点
- 精度高：两步法达到四阶精度，计算效率非常高
- 隐式方法：\(\beta_2 = 1/3 \neq 0\)，因此是隐式方法，需要迭代求解
- 基于Simpson求积公式：Milne方法实际上是用Simpson求积公式近似积分\(\int_{x_m}^{x_{m+2}} f(x, y(x)) dx\)得到的，这也是它四阶精度的来源
4. 稳定性问题
Milne方法虽然精度高，但存在一个严重的缺点：弱稳定性。它的绝对稳定区间非常小，甚至在某些情况下会出现寄生解增长的现象，因此不适合长时间积分。
为了克服Milne方法的稳定性问题，人们发展了Hamming方法，它也是四阶两步法，但稳定性好得多。
五、待定系数法的优势与局限性
1. 核心优势
- 通用性强：可以构造任意形式的线性多步法，而不仅仅是Adams类方法
- 灵活性高：可以通过调整自由参数，在精度、稳定性和计算量之间进行权衡
- 系统性强：有一套标准的流程来确定系数和分析精度
2. 局限性
- 阶数上限：根据Dahlquist第一障碍定理，k步显式线性多步法的最高精度阶数是k，k步隐式线性多步法的最高精度阶数是k+1（当k为奇数时）或k+2（当k为偶数时）
- 稳定性与精度的矛盾：通常精度越高的方法，稳定性越差，需要在两者之间进行平衡
六、线性多步法构造方法对比
| 构造方法 | 基本思想 | 适用范围 | 优点 | 缺点 |
|---|---|---|---|---|
| 插值积分法 | 用插值多项式近似被积函数，然后积分 | Adams类方法 | 直观易懂，物理意义明确 | 只能构造Adams类方法 |
| 待定系数法 | 通过Taylor展开比较系数确定参数 | 所有线性多步法 | 通用灵活，可构造任意形式的方法 | 计算量较大，物理意义不直观 |
七、核心知识点归纳总结
| 知识点 | 核心内容 | 关键公式/结论 |
|---|---|---|
| 一般形式 | \(\sum_{j=0}^k \alpha_j y_{m+j} = h \sum_{j=0}^k \beta_j f_{m+j}\) | \(\alpha_k=1\)，\(\beta_k=0\)为显式，否则为隐式 |
| 线性差分算子 | \(L[y(x);h] = \sum_{j=0}^k [\alpha_j y(x+jh) - h\beta_j y'(x+jh)]\) | |
| 阶数条件 | \(c_0=c_1=\dots=c_p=0, c_{p+1}\neq0\) | \(c_p = \frac{1}{p!}\sum j^p \alpha_j - \frac{1}{(p-1)!}\sum j^{p-1} \beta_j\) |
| 相容性条件 | \(c_0=0, c_1=0\) | \(\sum \alpha_j=0, \sum j\alpha_j = \sum \beta_j\) |
| 两步法一般公式 | 含一个自由参数\(\alpha\)的三阶方法 | 公式(1.4.36) |
| Milne方法 | \(\alpha=-1\)时的四阶两步法 | \(y_{m+2}=y_m+\frac{h}{3}(f_{m+2}+4f_{m+1}+f_m)\) 局部截断误差：\(-\frac{h^5}{90}y^{(5)}(\xi)\) |
学习建议
- 重点掌握待定系数法的步骤：这是构造线性多步法的通用方法，必须熟练掌握从算子定义到阶数条件推导的完整流程。
- 理解线性差分算子的作用：算子是将微分方程离散化的数学工具，理解它有助于掌握数值方法的本质。
- 记住Milne方法的公式和特点：Milne方法是最著名的线性多步法之一，虽然稳定性有缺陷，但它的精度和效率非常高。
- 了解阶数上限定理：知道线性多步法的精度不是无限提高的，存在理论上的上限。
posted on 2026-05-09 21:39 Indian_Mysore 阅读(9) 评论(0) 收藏 举报
浙公网安备 33010602011771号