符号说明
- - 微分方程的精确解（理论解）函数
- - 在点 处的数值解（近似值）
- - 在点 处的数值解（近似值）
- - 函数 y 的一阶导数，等于
- - 函数 y 的二阶导数，用于泰勒展开分析误差
- - 微分方程右端的已知函数，表示导数 的计算方式
- - 第 n 个离散点的横坐标，数值解的计算点
- - 第 n+1 个离散点的横坐标，下一个计算点
- - 步长，即 ，影响计算精度和效率
- - 介于 和 之间的某个值（用于泰勒展开余项）
- - 表示误差阶数为一阶，如欧拉法的全局截断误差
- - 表示误差阶数为二阶，如欧拉法的局部截断误差和梯形法的全局截断误差
常微分方程一般形式
常微分方程的一般形式可以表示为：
其中：
- y 是未知函数
- x 是自变量
- 是已知函数
给定初值条件：
这种形式的常微分方程称为初值问题或柯西问题。
计算目标
- 数学家：求出y的表达式
- 工程师：给定x，y大概是多少？
常微分方程数值解法的基本逻辑
常微分方程数值解法的基本思路是：已知自变量处的函数值，利用微分方程和数值计算方法求解下一个点处的函数值。
其中，步长是一个关键参数，定义为：
步长的选择需要权衡计算精度和计算效率：
- 较小的步长通常可以获得更高的计算精度，但会增加计算量
- 较大的步长可以减少计算量，但可能导致累积误差增大
通过不断重复这个推进过程，并合理选择步长，我们可以逐步获得微分方程在一系列离散点上的数值解。不同的数值方法主要区别在于如何由推导出的具体计算方案。
（向前）欧拉法
欧拉法是最简单的常微分方程数值解法之一。我们可以从泰勒级数展开的角度来推导欧拉法。
设在点处有连续导数，则在附近的泰勒展开式为：
其中是介于和之间的某个值。
令为步长，并截断二阶及以上的余项，得到：
由于，因此可以得到欧拉法的递推公式：
这就是欧拉法的基本形式。它是一种一阶精度的显式方法，局部截断误差为，全局截断误差为。
欧拉法的几何意义
欧拉法的几何意义可以理解为：在每一步计算中，都是用该点处的切线方向来近似函数在该区间的变化。具体来说：
- 在点处，斜率确定了一条切线
- 沿着这条切线方向前进步长，得到下一个点
- 重复这个过程，不断用局部的切线来逼近真实的解曲线
这种方法的误差主要来源于：实际的解曲线是弯曲的，而我们用分段直线（切线）来逼近它。这也解释了为什么步长越小，通常可以得到更准确的结果——因为用更短的切线段来逼近弯曲的曲线，累积的误差会更小。
欧拉法的数值积分推导
欧拉法也可以从数值积分的角度来理解和推导。考虑微分方程：
对等式两边在区间上积分：
左边积分后得到：
右边积分可以用左矩形法（最简单的数值积分方法）来近似：
将上述结果代入，得到：
这与我们之前从泰勒展开得到的欧拉法公式完全一致。这种推导方式突出了欧拉法的另一个解释：它实际上是在用最简单的数值积分方法（矩形法）来计算导数的积分。
这也解释了为什么欧拉法的精度较低——因为矩形法是最基本的数值积分方法，其精度相对较低。
欧拉法的误差与精度
从泰勒展开的角度来看，在点处展开到二阶项：
欧拉法只保留到一阶项：
因此，每一步的局部截断误差（Local Truncation Error, LTE）为泰勒展开中被截断的第一项：
在区间上，总共需要步，每一步都会产生的局部误差。这些误差会累积，导致全局截断误差（Global Truncation Error, GTE）为：
因此，欧拉法的全局截断误差为，即为一阶精度方法。这表明当步长减小为原来的一半时，理论上误差也会减小为原来的一半。
例1：
例2：试用欧拉法计算下列初值问题：
取步长 ，并与精确值 比较。
解：由欧拉法得：
,
计算结果表格：
精确值 | 误差 | |||
0 | 0 | 0 | 0 | 0 |
1 | 0.5 | 0.500000 | 0.433333 | 0.066667 |
2 | 1.0 | 0.800000 | 0.666667 | 0.133333 |
3 | 1.5 | 0.900000 | 0.807692 | 0.092308 |
4 | 2.0 | 0.984615 | 0.933333 | 0.051282 |
每一步算的都是不同的值；因此误差并不会随着迭代收敛；类似随机出现
后退欧拉法（隐式方法）
后退欧拉法（Backward Euler Method）是欧拉法的一个变体，也称为隐式欧拉法。与前向欧拉法不同，它使用终点处的函数值来计算数值解。其基本公式为：
注意到公式右边包含了未知量，这使得后退欧拉法成为一个隐式方法。
- 通常需要通过迭代方法（如牛顿法）来求解每一步的。
- 考试中，通过化简一般是可以得到关于……的显示表达式
后退欧拉法的积分推导
类似于前向欧拉法，我们也可以从积分的角度推导后退欧拉法。考虑微分方程：
在区间上积分：
左边积分得到：
与前向欧拉法不同，后退欧拉法使用右矩形法来近似右边的积分：
代入得到：
这就是后退欧拉法的基本公式。与前向欧拉法相比，它使用了区间右端点的函数值，这使得方程变成了一个隐式方程。
梯形法（隐式方法）
梯形法是一种数值求解常微分方程的方法，它结合了前向欧拉法和后退欧拉法的思想，使用了积分区间两端点的函数值的平均值来近似积分。其基本公式为：
梯形法的积分推导
这个公式可以从积分的角度来理解。考虑微分方程：
在区间上积分：
左边积分得到：
右边的积分可以用梯形法来近似：
这是一个隐式方法，因为公式右边包含了未知量。与欧拉法相比，梯形法具有更高的精度，其局部截断误差为，全局截断误差为。
梯形法的优点：
- 更高的精度：是二阶方法
- 更好的稳定性：比欧拉法具有更好的数值稳定性
- 对称性：同时考虑了区间两端点的函数值
梯形法的缺点：
- 计算复杂：每一步都需要求解非线性方程
- 计算成本：比显式方法（如欧拉法）需要更多的计算
例3：试用梯形法计算下列初值问题：
取步长 ，并与精确值 比较。
解：由梯形法得：
化简得：
计算结果表格：
精确值 | 误差 | |||
0 | 0 | 0 | 0 | 0 |
1 | 0.5 | 0.416667 | 0.433333 | 0.016666 |
2 | 1.0 | 0.666667 | 0.666667 | 0 |
3 | 1.5 | 0.812500 | 0.807692 | 0.004808 |
4 | 2.0 | 0.937500 | 0.933333 | 0.004167 |
改进欧拉公式（预报-校正格式）
改进欧拉法（Modified Euler Method）是一种预报-校正格式的数值方法，它通过两步计算来提高精度：预报步和校正步。
基本思想
- 预报步：使用前向欧拉法计算一个初步预测值
- 校正步：使用预测值来计算最终的数值解
其中：
- 是预报步得到的预测值
- 是校正步得到的最终数值解
- 最终表达式
将预测值代入校正步公式：
这就是改进欧拉法的最终表达式，它是一个显式公式，可以直接计算 。
改进欧拉法的特点
- 精度：局部截断误差为，全局截断误差为
- 稳定性：比简单欧拉法具有更好的稳定性
- 计算量：每步需要计算两次函数值
- 实现简单：不需要求解非线性方程（与隐式方法相比）
改进欧拉法可以看作是介于显式欧拉法和梯形法之间的一种方法。它保持了显式方法的简单性，同时提高了计算精度，精度在欧拉法与梯形法之间。
龙格-库塔法(R-K法)
建立高精度的单步递推公式
龙格-库塔法的基本思想
龙格-库塔法的基本思想来源于微分中值定理。根据微分中值定理，在区间上存在一点，使得：
问题的关键在于如何找到这个合适的中间点，或者说如何得到亦即的一个良好近似。龙格-库塔法就是通过在区间内选取若干个精心设计的点，对其斜率进行加权平均来近似这个导数值。
本质上，目前所有方法都是在逼近微分中值定理的那个斜率
四种方法的逻辑
从微分中值定理的角度来看，这些数值方法的本质区别在于如何选取斜率来近似区间内的平均变化率：
- 向前欧拉法：直接使用区间左端点的斜率 ，最简单但精度较低
- 向后欧拉法：使用区间右端点的斜率 ，需要迭代求解隐式方程
- 梯形法：取区间两端点斜率的算术平均值 ，是一种隐式方法
- 改进欧拉法：先用向前欧拉法预测右端点值，再用这个预测值计算右端点斜率，最后取两个斜率的平均值，避免了求解隐式方程
- 龙格库塔法：在区间内选取多个精心设计的点计算斜率，通过加权平均得到更好的近似，是这些方法中最精确的
从计算复杂度来看，这些方法呈现出一个渐进的过程：向前欧拉法最简单但精度最低，龙格库塔法最复杂但精度最高。每种方法都可以看作是对前一种方法的改进，通过更复杂的斜率选取方式来提高精度。
二阶龙格库塔
二阶龙格库塔法是在泰勒展开的基础上构造的一种数值解法。我们从以下步骤推导：
1. 基本形式
假设二阶龙格库塔法具有如下形式：
其中、、是待定系数，满足：
2. 泰勒展开
为了详细理解二阶龙格库塔法的精度，我们将真解在点处展开至二阶项：
同时，将数值解展开：
其中：
将和代入得到：
与真解对比，通过配系数得到：
3. 最终形式
一组常用的系数选择是：
这样得到的二阶龙格库塔法（中点法）为：
这种方法的局部截断误差为，全局截断误差为。
另一种常用的参数选择是：
这样得到的是改进欧拉法：
这与我们之前介绍的改进欧拉法完全一致，说明改进欧拉法实际上是二阶龙格库塔法的一个特例。
我们实际上只能控制α，也就是只能控制第二个采样点的位置
4. 各参数的意义
在二阶龙格库塔法中，这三个参数具有以下具体含义：
- ：控制第二个采样点的位置，表示在区间内向前推进的比例。
- 当 时，在区间中点取样（中点法）
- 当 时，在区间终点取样（改进欧拉法）
- 和 ：权重系数，用于对两个斜率进行加权平均
- 满足 ，确保加权和为1
- 表示起点斜率的权重
- 表示第二个采样点斜率的权重
这些参数的选择直接影响方法的稳定性和精度。例如，中点法（）和改进欧拉法（）都是通过精心选择这些参数得到的特殊情况。
三阶龙格-库塔
三阶龙格-库塔法的一般形式为：
约束条件：
为了保证三阶龙格-库塔法具有三阶精度，参数需要满足以下约束条件：
以上约束条件各自有其特定的作用：
- ：确保权重系数之和为1，这是加权平均的基本要求
- ：保证第三个采样点的位置与其计算中使用的系数一致性
- ：确保方法具有二阶精度，与泰勒展开二阶项系数匹配
- ：保证三阶项系数的准确性，是实现三阶精度的必要条件
- ：与三阶导数项相关，用于消除三阶截断误差
这些条件共同作用，确保了方法具有三阶精度。
可以看到，约束条件是很多的；虽然参数有无数种组合，但是参数的选择并不是任意的
最终形式
一种常用的三阶龙格-库塔法参数选择是：
代入上述参数后得到具体计算公式：
相当于选取了：起点、中点、终点三个点，且中点的权重最大
这种方法的局部截断误差为，全局截断误差为。
四阶龙格-库塔
四阶龙格-库塔法的一般形式为：
约束条件：
为了保证四阶龙格-库塔法具有四阶精度，参数需要满足以下约束条件：
这些约束条件各自有其特定的作用：
- ：确保权重系数之和为1，这是加权平均的基本要求
- ：保证第三个采样点的位置与其计算中使用的系数一致性
- ：保证第四个采样点的位置与其计算中使用的系数一致性
- ：确保方法具有二阶精度
- ：确保方法具有三阶精度
- ：确保方法具有四阶精度
- ：与高阶导数项相关，用于消除更高阶的截断误差
这些条件共同作用，确保了方法具有四阶精度。
最终形式
通过与泰勒展开对比可以得到一系列条件方程。最常用的参数选择是：
代入这组参数后得到经典的四阶龙格-库塔法：
这种方法的局部截断误差为，全局截断误差为。它在精度和计算效率之间取得了很好的平衡，是求解常微分方程最广泛使用的数值方法之一。
例：试用数值方法计算下列初值问题：
分别应用欧拉法、后退欧拉法、梯形法和龙格-库塔法。
多步法
单步法与多步法
计算方法基本观念：用更多、更新的数据得到的数值解法精度更高
单步法与多步法的主要区别：
- 信息利用：
- 单步法：仅使用当前点的信息计算下一点
- 多步法：利用多个历史点的信息计算下一点
- 计算特点：
- 单步法：每步计算独立，需多次计算函数值（如龙格-库塔法）
- 多步法：可重复利用历史数据，减少函数求值次数
- 启动过程：
- 单步法：只需初值即可直接启动
- 多步法：需要其他方法（通常是单步法）提供初始几个点的值
即使是RK法，在一个步长内精心选取了若干个点，但是实际上还是只用到了上一个数据点的信息，所以还是单步法
多步法的基本思路
多步法是一类利用多个已知点的信息来计算下一个点的数值解法。与单步法（如龙格-库塔法）相比，多步法可以更有效地利用已计算的信息。
基本特点
- 利用已知的多个点的信息进行计算
- 可以减少函数求值次数，提高计算效率
- 需要其他方法（如龙格-库塔法）来获取初始值
一般形式
可以看作函数值和切线斜率的线性组合
其中：
- 表示使用的历史点数（影响方法的阶数）
- 为待定系数，一共有个（通过与泰勒展开对比确定）
- 表示在点 处的函数值
使用历史点数多→待定系数多→解出待定系数需要的方程多→泰勒展开阶数高
约束条件
为了确定多步法中的系数，我们需要将理论解和近似解在点 处展开成泰勒级数并进行比较。
- 基本数值解泰勒级数表达式：
- 解析解泰勒级数展开式：
- 比较两式，对应系数相等得到约束条件：
（3）和（4）一起，一共是个方程；m取决于泰勒展开到哪一阶，理论上可以无限大
多步法的底层逻辑是：用更多的历史的点就能算得更准；Adam法关注b，把a全部置零；Gear法关注a，把b置零
Adam法
从约束条件可以看出：
- 我们一共有个方程（一个基本约束方程和个高阶项约束方程）
- 但未知系数的个数为个（个和个系数）
虽然理论上泰勒级数可以无穷展开，m可以取无穷大，但是这样计算压力很大（方程阶数提高）
当时，方程组有无穷多组解。这给了我们选择系数的自由度，可以根据不同的需求选择不同的系数组合。Adam方法就是一种特殊的系数选择方案。
相当于方程数目不够，我们额外给定一些条件
- 特殊系数选择(p个额外条件，所有Adam方法都有这个操作)
- 简化后的多步法公式：
显式Adam法（Adam's-Bashforth法）：
也就意味着不会用到n+1的条件
结合式6可知，这意味着令= 0；由（4）式得到约束方程：
- 求解条件：
至此，我们额外给定了个条件
还剩下：
个未知数，个方程，要求
式7也可写成矩阵形式：
也可写为更好记忆的形式：
隐式Adam法（Adam-moulton法）
此时
p+2个未知数，m个方程，所以应有：
Adams隐式方法的矩阵方程形式：
也可写为：
式（12）和式（11）本质上完全相同，几乎只用背一个（1~0~ -1~-2~……-p）
Gear法
Gear法：除外所有的都为零，因为，所以所有的Gear法都为隐式法。
k阶Gear公式的推导为：
令p=k-1和，得k阶Gear法：
与Adam法类似，式(11)中的k+1个系数可由下式得出：
数值稳定性
数值稳定性的数学定义
对于常微分方程初值问题：
如果在求解过程中出现的任何扰动（包括舍入误差和截断误差）在后续的步骤中能够逐步衰减，则称这种方法是稳定的
测试方程法
探究步长和稳定性的关系
测试方程法是一种评估数值方法稳定性的重要手段。其基本思想是选取一个简单的线性测试方程来分析数值解法的稳定性特征。
其中λ是一个复数。这个方程的精确解为：
选择这个测试方程的原因：
- 简单性：方程形式简单，便于分析
- 代表性：许多实际问题在线性化后都可以归结为这种形式
- 特征明显：当Re(λ) < 0时，精确解是衰减的；当Re(λ) > 0时，精确解是增长的
稳定性分析步骤：
- 将数值方法应用于测试方程
- 得到差分方程
- 分析差分方程的特征方程
- 确定稳定性条件
通过测试方程法，我们可以得到数值方法的稳定域，即保证数值解稳定的所有复数λh的集合。
对测试方程应用欧拉法
- 将前向欧拉法应用于测试方程：
- 差分方程形式：
- 稳定性条件：
- 递推关系分析：由于 ，这是一个等比数列。如果要使数值解稳定（不发散），必须保证公比的绝对值不大于1，即
- 误差分析：假设在某一步引入扰动 ，则下一步的扰动为 。要使误差不放大，必须满足
- 与精确解比较：测试方程的精确解为 。当 时精确解是衰减的，数值解也应该表现出相同的衰减特性，这要求
我们可以从以下几个角度理解这个稳定性条件：
- 几何解释：
- 在复平面上，必须落在以点为中心，半径为1的圆内部及边界
- 这表明欧拉法是条件稳定的
- 意味着λ越大，步长h必须选择得越小
对测试方程应用后退欧拉法
- 将后退欧拉法应用于测试方程：
- 整理得到差分方程形式：
- 稳定性条件：
- 几何解释：
- 在复平面上，hλ必须位于以(1,0)为圆心的单位圆外
- 这表明后退欧拉法是无条件稳定的
三阶Gear法的绝对稳定域
三阶Adam-moulton的绝对稳定域
Adam-Bashforth法稳定性
阶数越高，对稳定性要求越苛刻
隐式方法稳定性优于显式方法
、
