——————————————————————————————————————————————
数值分析笔记
——————————————————————————————————————————————
Tsui Dik Sang
2024 年 2 月 18 日——2025 年 5 月 15 日
数值分析笔记 Tsui Dik Sang
2
第一章 写在笔记之前
虽然老师说数值分析不算是数学课,然而我觉得如果能尽量理解清楚每一个方法的数学原理,对于打通整本书的逻辑,并灵
活应用这些方法都是有益的,这本笔记也正是记录了笔者尝试用自己的语言来将课本一些较为抽象但是最终还是消化理解了逻辑
的定理描述出来。
当然,并非所有逻辑都搞清了,仍然有大量的内容没有彻底弄清楚,比如课程规划未讲授的第三章和第八章,以及共轭梯度
法、CG 等,在已学过的内容中也有一些因时间精力受限而只是囫囵吞枣记住结论没有理解推导的定理,对于这类“黑箱”,我尝
试尽可能的将其封装——我不需要知道里面是什么,如果用一个芯片去比喻的话,我希望看到他的输入输出引脚都暴露出来,而
内部结构完全黑箱化,或者可能会暴露一些电路,我只需要知道他需要我输入什么,他能给我输出什么,至于他可以通过某个函
数变换和公式实现这样的效果,我不管,正如我不需要知道这些电路是如何想到如何搭建设计的,但是我至少还是需要了解一些
修电路的知识,以便能在芯片发生不大的可逆故障的时候可以手动调整修复,或者根据自己的使用环境做适当的优化布局。
最后,考试不是终点,学数值分析更加依赖于编程和查表,不一定每一个方法都要清晰记住,但是需要有属于自己的检索系
统,并且电子的最优,以便能在需要用到的时候快速锁定最适的方法,并且有现成的实现代码模版,那就相当于一个芯片,焊在
主板上就能实现一个强大的扩展功能了,我将自己上课和作业的一些代码整理到了仓库 MycodeforNumericalAnalysis,,尽量做
了函数的封装,当然,也不算完善,只有有时间会补上,也欢迎有兴趣的一起来贡献!
笔记封面的 Lorenz 吸引子是我在课堂报告中研究的一个现象。我发现这个动力系统非常有趣。记得当时,我对着这个图
像反复调整参数,观察着形成的非常具有“立体感”的图形,这种立体感比正方体球等立方体更加强烈,是难以用二维的
视角去想象的,因为我似乎无法预测我转动一个角度后呈现给我的二维剖面是什么,不知不觉就沉浸其中很长时间——这
大概就是数值计算之美的一个缩影吧。
Tsui Dik Sang
—2024.2.18
3
数值分析笔记 第一章 写在笔记之前 Tsui Dik Sang
4
目录
第一章 写在笔记之前 3
1.1 误差 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
1.1.1 误差的分类 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
1.1.1.1 误差来源 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
1.1.1.2 绝对误差 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
1.1.1.3 相对误差 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12
1.1.2 有效数 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12
1.1.3 误差的估计 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13
1.1.3.1 基本误差估计 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13
1.1.3.2 误差的传播 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13
1.2 误差定性分析 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14
1.2.1 算法的稳定性 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14
1.2.1.1 算法 A: 正向递推公式法 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14
1.2.1.2 算法 B: 反向递推法 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
1.2.1.3 稳定的算法 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
1.2.2 病态问题判断 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
1.2.3 避免误差伤害 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16
1.2.4 算法设计 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16
1.2.5 加权平均的松弛技术 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16
第二章 插值法 17
2.1 插值问题与多项式插值简介 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17
2.1.1 插值问题 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17
2.1.2 多项式插值 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17
2.1.2.1 多项式插值的唯一性 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
2.2 拉格朗日插值法 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
2.2.1 一次插值 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
2.2.2 Lagrange 插值公式 (多次插值) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
2.2.3 插值余项与误差估计 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19
2.2.3.1 用零点式 ω 表示的插值函数 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19
2.2.3.2 定理阐释 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19
2.2.3.3 定理证明 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20
2.2.3.4 余项的估计 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20
2.2.3.5 余项定理的应用 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21
2.3 Newton 插值法 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21
2.3.1 插值多项式的逐次生成以及均差的引入 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21
5
数值分析笔记 目录 Tsui Dik Sang
2.3.1.1 一次差值 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21
2.3.1.2 二次差值 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21
2.3.1.3 均差的引入与推广到多次 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22
2.3.1.4 Newton 均差插值多项式 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23
2.3.2 差分形式的 Newton 插值多项式 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23
2.3.2.1 差分算子 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23
2.3.2.2 差分算子的运算 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24
2.3.2.3 Newton 前插公式 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25
2.4 Hermite 插值 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25
2.4.1 函数导数的均差表示 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25
2.4.1.1 重节点均差 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25
2.4.2 Taylor 插值多项式 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26
2.4.2.1 多项式得到 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26
2.4.2.2 余项误差 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26
2.4.2.3 一般规律总结 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26
2.4.3 Hermite 插值插值:利用重节点均差 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27
2.4.3.1 待定系数法 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28
2.4.3.2 重节点均差法 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28
2.4.4 两个典型的 Hermite 插值 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28
2.4.4.1 P (x
i
) = f(x
i
)(i = 0, 1, 2) & P
′
(x
1
) = f
′
(x
1
) . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28
2.4.4.2 两点三次多项式插值 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29
2.5 分段低次插值 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29
2.5.1 高次插值的病态性质 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29
2.5.2 分段线性插值 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30
2.5.3 分段三次插值 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30
2.6 三次样条插值 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31
2.6.1 条件的确定 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31
2.6.1.1 罗列已知条件 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31
2.6.1.2 对比目标所需条件数 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31
2.6.1.3 尝试构造缺失的条件 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31
2.6.2 样条插值函数的建立 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32
2.6.2.1 第一类边界条件:已知边界一阶导 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32
2.6.2.2 第二类边界条件:已知边界二阶导 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32
2.6.2.3 第三种边界条件:周期 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33
第三章 线性方程组解法 35
3.1 引言:一些基本的矩阵概念 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35
3.1.1 非奇异矩阵 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35
3.1.2 特征值与谱半径 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35
3.1.3 Jordan 标准型 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35
3.1.4 范数 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36
3.1.4.1 向量的范数 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36
3.1.5 矩阵范数 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36
3.1.6 向量、矩阵序列的极限 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37
3.1.6.1 基础定义和结论 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37
3.1.6.2 极限与谱半径 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37
6
数值分析笔记 目录 Tsui Dik Sang
3.1.7 对角占优矩阵等矩阵 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38
3.1.7.1 定义 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38
3.1.7.2 定理 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38
3.2 直接解法 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39
3.2.1 高斯消元法 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39
3.2.1.1 列主元高斯消元法 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39
3.2.2 三角分解法 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40
3.2.2.1 LU 分解 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40
3.2.2.2 平方根法 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40
3.2.2.3 追赶法 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40
3.2.3 误差分析 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41
3.2.3.1 b 的变换对解的影响 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41
3.2.3.2 A 的变换对解的影响 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41
3.2.4 条件数 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41
3.3 迭代法 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42
3.3.1 迭代法的思想 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42
3.3.1.1 算法 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42
3.3.1.2 收敛性 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42
3.3.2 Jacobi 迭代法 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42
3.3.2.1 基本定义 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42
3.3.2.2 收敛性 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43
3.3.3 Gauss Seidel 迭代法 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43
3.3.3.1 基本定义 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43
3.3.3.2 收敛性 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44
3.4 超松弛迭代法 (SOR) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45
3.4.1 方法 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45
3.4.2 收敛性 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45
第四章 非线性方程数值解法 47
4.1
基础定义和定理
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47
4.2 二分法 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47
4.2.1 近似程度 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47
4.2.2 优点和缺点 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48
4.2.2.1 优点 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48
4.2.2.2 缺点 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48
4.3 不动点法 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48
4.3.1 算法简述 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48
4.3.2 性质 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49
4.3.2.1 唯一性 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49
4.3.2.2 收敛性 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50
4.3.2.3 误差估计公式 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50
4.3.2.4 迭代函数在解处的导数 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50
4.3.3 迭代法的收敛阶 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50
4.3.3.1 收敛阶的定义和一些判定 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50
4.3.3.2 两种算法的收敛阶讨论 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51
4.3.4 优化的迭代法 * . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51
7
数值分析笔记 目录 Tsui Dik Sang
4.3.4.1 线性收敛序列的 Aitken 加速法 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51
4.3.4.2 Steffensen 迭代法 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52
4.4 Newton 迭代法 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52
4.4.1 简化 Newton 法:平行弦方法 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52
4.4.2 Newton 下山法 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53
4.4.3 重根情况 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53
4.4.3.1 重根情形分析 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53
4.4.3.2 平方收敛的重根算法 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53
4.4.4 类差分插值法 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54
4.4.4.1 弦截法 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54
4.4.4.2 抛物线法 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54
4.5 非线性方程组的数值解法 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54
4.5.1 非线性方程组的矩阵表示 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54
4.5.2 不动点迭代法 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55
4.5.2.1 基本使用 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55
4.5.2.2 误差分析 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55
4.5.2.3 收敛性 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55
4.5.3 Newton 迭代法 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55
4.5.3.1 基本使用 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55
4.5.3.2 收敛性 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56
第五章 数值积分与数值微分 57
5.1 数值积分的基础概念 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57
5.1.1 思想与引入 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57
5.1.1.1 两点积分近似 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57
5.1.1.2 多点扩展 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57
5.1.2 代数精度 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57
5.1.3 求积公式的余项 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58
5.1.4 收敛性与稳定性 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58
5.1.4.1
收敛性
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58
5.1.4.2 稳定性 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58
5.2 插值型求积 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59
5.2.0.1 求积公式 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59
5.2.0.2 余项公式 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59
5.3 Newton-Cotes 公式 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59
5.3.1 引入 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59
5.3.2 余项 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60
5.3.2.1 奇数阶 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60
5.3.2.2 偶数 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60
5.3.3 Simpson 公式 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60
5.3.3.1 定义 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60
5.3.3.2 余项 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61
5.4 复合求积公式 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61
5.4.1 复合梯形公式 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61
5.4.1.1 余项 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61
5.4.2 复合 Simpson 公式 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61
8
数值分析笔记 目录 Tsui Dik Sang
5.4.2.1 余项 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62
5.5 Romberg 求积公式 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62
5.5.1 主线:梯形公式的递推化 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62
5.5.2 支线:外推加速算法 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63
5.6 自适应算法 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64
5.7 Gauss 求积公式 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64
5.7.1 二点插值 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64
5.7.2 权函数 ρ(x) 的引入 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65
5.7.3 Gauss-Legendre 公式 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65
5.8 多重积分 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66
5.9 数值微分 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66
5.9.1 中点方法 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66
5.9.2 插值型求导公式 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66
第六章 常微分方程初值问题数值解法 67
6.1 引入 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67
6.1.1 单步法 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67
6.1.2 Euler 法 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67
6.1.2.1 方法推导 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67
6.1.2.2 后退的 Euler 法 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68
6.1.2.3 误差 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68
6.1.3 梯形方法 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68
6.1.4 改进的 Euler 法 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68
6.1.5 单步法总结及误差分析 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69
6.2 Runge-Kutta 方法 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70
6.2.1 一般形式 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70
6.2.2 二阶显式 Runge-Kutta 法 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71
6.2.2.1 基本定义 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71
6.2.2.2 常数的确定 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71
6.2.3
三阶、四阶显式
Runge-Kutta
法
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72
6.2.3.1 三阶显式 Runge-Kutta 法 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72
6.2.3.2 四阶显式 Runge-Kutta 法 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72
6.2.4 变步长 Runge-Kutta 法 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73
6.3 单步法误差分析 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73
6.3.1 收敛性 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73
6.3.1.1 基础定义 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73
6.3.2 相容性 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74
6.3.3 稳定性 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75
6.3.3.1 模型方程 y
′
= λy 的建立 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75
6.3.3.2 Euler 法的稳定性 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75
6.3.3.3 二阶 Runge-Kutta 法的稳定性 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76
6.3.3.4 更高阶的 Runge-Kutta 法的稳定性 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76
6.3.3.5 隐式 Euler 法的稳定性 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76
6.3.3.6 梯形法的稳定性 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76
6.4 线性多步法 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77
6.4.1 一般公式与基础概念 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77
9
数值分析笔记 目录 Tsui Dik Sang
6.4.2 Adams 显式与隐式多步法 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78
6.4.3 更多经典线性多步法 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79
6.4.3.1 Milne 法 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79
6.4.3.2 Simpson 法 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79
6.4.3.3 Hamming 法 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80
6.4.3.4 PECE 预测-校正法 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80
6.4.3.5 PMECME 预测-校正法 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80
6.4.3.6 总结 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81
6.4.4 收敛性与稳定性 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81
6.4.4.1 相容性 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81
6.4.5 一阶方程组域刚性方程组 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81
6.4.5.1 一阶方程组 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81
6.4.5.2 化高阶方程为一阶方程组 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82
6.4.5.3 刚性方程组 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82
第七章 笔记初稿完成 85
10
概论:误差与算法
的无他,只是给出了一些基本的定义之类的东西咯,也是对数值计算这一门课的一个范围界定。
1.1 误差
1.1.1 误差的分类
1.1.1.1 误差来源
1
定义 1.1.1 (观测误差). 观测得到的物理类所存在的误差
定义 1.1.2 (截断误差/方法误差). 用数值方法求解近似解时,近似解与精确解的误差就是截断误差/方法误差
例如用泰勒多项式逼近可微函数
R
n
(x) = f(x) − P
n
(x) =
f(n + 1)(ξ)
(n + 1)!
x
n+1
(1.1)
就是截断误差。
这本书所讨论的误差基本都是截断误差
1.1.1.2 绝对误差
定义 1.1.3 (绝对误差). 若 x
∗
为 x 的近似值,那么绝对误差为
e(x
∗
) = x
∗
− x (1.2)
由于实际的真值也是未知的,只能根据理论计算估计,所以绝对误差也不能描述清楚,只知道一个限度
定义 1.1.4 (绝对误差限 = 误差). 可以估计出
|e(x
∗
)| = |x
∗
− x| ⩽ ε (1.3)
记做 x = x
∗
± ε
为了消除量纲对刻画近似准确度的估计,引入了相对误差
1
误差来源是 [2] 这本书上补充的
11
数值分析笔记 目录 Tsui Dik Sang
定义 1.1.5 (相对误差). 定义相对误差为
e + r(x
∗
) =
x
∗
− x
x
∗
(1.4)
1.1.1.3 相对误差
同样也可以定义相对误差限
定义 1.1.6 (相对误差限).
ε
r
(x
∗
) =
ε(x
∗
)
|x
∗
|
(1.5)
1.1.2 有效数
引入有效数字以及有效数, 先直接给书中的严谨定义
2
定义 1.1.7 (有效数字). 设近似值可以表示成下面的形式
x
∗
=
p
X
i=1
10
m−i
a
i
= ±10
m
× 0.a
1
a
2
···a
n
···a
p
, a
1
̸= 0 (1.7)
若使得不等式
|
x
∗
−
x
|
⩽
10
m−n
2
(1.8)
成立的最大整数为 n,则称近似值 x
∗
具有 n 位有效数字,分别是 a
1
, a
2
, ··· , a
n
定义 1.1.8 (有效数). 特别的,如果 n = p, 那么 x
∗
被称为有效数
上面从公式上的定义可能容易数晕位数,那么我们再给出书 [2] 中给的文字定义
定义 1.1.9 (有效数字的文字定义). 若近似值 x
∗
的误差限是某一位的半个单位,该位到 x
∗
的第一位非零数字共有 n 位,那
么称 x
∗
有 n 位有效数字,
上面的这个定义实际上是从误差的角度定义了有效数字,即
有效数字的本质实际是一种误差
其实我们之前接触过有效数字的通俗定义
3
定义 1.1.10 (有效数字的通俗定义). 假设 x
∗
是 x 的近似值,如果从第一个非零数字开始的第 n 位对 x
∗
和 x 进行四舍五入
2
这里给的定义是 [1] 中的说法,[2] 中用
的形式来表示误差,这里需要尤其注意小数点的位置,才能更好的理解式1.8,书 [2] 中是形如 10
m
×a
n
×10
−(n−1)
的对应形式,那么对应的误差不等式就应该是
|x
∗
− x| ⩽
10
m−n+1
2
(1.6)
反正注意判断
3
源于csdn 某博客启发
12
数值分析笔记 目录 Tsui Dik Sang
结果完全一致,那么称 x
∗
有 n 位的有效数字
将定义1.1.10与定义1.1.7对比着看就很清晰的看到式1.8的直观意义。
在 [2] 中引入了相对误差限的定理
4
定理 1.1.1 (相对误差限与有效数字). 若 x
∗
有 n 位的有效数字,那么
ε
∗
r
⩽
1
2a
1
× 10
−(n−1)
(1.11)
以及其逆定理
5
推论 1.1.2 (由相对误差限判断有效数字). 若 x
∗
的相对误差限满足
ε
∗
r
⩽
1
2(a
1
+ 1)
× 10
−(n−1)
(1.13)
那么 x
∗
至少有 n 位有效数字
1.1.3 误差的估计
1.1.3.1 基本误差估计
实际上指的是基本运算的误差传播
6
,只有几个公式
定理 1.1.3 (四则运算的误差传递公式).
ε
(
x
∗
1
± x
∗
2
) ⩽ ε(x
∗
1
) + ε(x
∗
2
) (1.14)
ε(x
∗
1
x
∗
2
) ⩽ |x
∗
1
|ε(x
∗
2
) + |x
∗
2
|ε(x
∗
1
) (1.15)
ε
x
∗
1
x
∗
2
⩽
|x
∗
1
|ε(x
∗
2
) + |x
∗
2
|ε(x
∗
1
)
|x
∗
2
|
2
, x
∗
2
̸= 0 (1.16)
1.1.3.2 误差的传播
定理 1.1.4 (多元函数误差的传播). 假设有 n 个自变量 x
1
, x
2
, ··· , x
n
,通过函数 f 构成因变量 y = f (x
1
, x
2
, ··· , x
n
),那么
因变量的误差与自变量的误差关系为
e(y
∗
) = y
∗
− y ≈
n
X
i=1
f
i
(x
∗
1
, x
∗
2
, ··· , x
∗
n
)e(x
∗
i
) (1.17)
4
注意到下面不等式
a
1
× 10
m
⩽ |x
∗
| ⩽ (a
1
+ 1) × 10
m
(1.9)
当 x
∗
有 n 位的有效数字时, 结合1.9与1.8, 即可推出
ε
∗
r
=
ε
∗
|x
∗
|
⩽
0.5 × 10
m−n+1
a
1
× 10
m
=
1
2a
1
× 10
−n+1
(1.10)
5
同理,用1.9来证
|x − x
∗
| ⩽ ε
∗
= |x
∗
|ε
∗
r
< (a
1
+ 1) × 10
m
×
1
2(a
1
+ 1)
× 10
−n+1
(1.12)
满足1.8, 因此 x
∗
至少有 n 位有效数字。
6
更一般的会在“误差的传播”提到
13
数值分析笔记 目录 Tsui Dik Sang
其中 f
i
=
∂f
∂x
i
这个定义的证明不难,用的是微分中值定理,因此中间实际上用的是约等于好,并不是全等。假设除了 x
i
之外都没有误差,推
出一元的公式,同理推出 n 个一元结果,再求和即可。
同样有相对误差的递推公式
推论 1.1.5 (相对误差的递推公式).
e
r
(y
∗
) =
n
X
i=1
f
i
(x
∗
1
, x
∗
2
, ··· , x
∗
n
)
x
∗
i
y
∗
e
r
(x
∗
i
) (1.18)
1.2 误差定性分析
7 8
1.2.1 算法的稳定性
其实也就是讨论初始数据的误差在误差传播中对最后结果的影响效果,下面用一个例子来说明这个影响到底有多大。
计算 I
n
= e
−1
ˆ
1
0
x
n
e
x
dx,并估计误差
1.2.1.1 算法 A: 正向递推公式法
由分部积分可以得到
I
n
= 1 − nI
n−1
, n = 1 , 2, ···
I
0
= e
−1
ˆ
1
0
e
x
dx = 1 − e
−1
(1.19)
然后我们用泰勒公式
9
来确定初值, 就可以构造算法 具体的结果不放了,但是在取 k=7、n=8 的时候出现了负值,这是一个很大
算法 1 算法 A
使用1.20确定初值
˜
I
0
for all n = 1, 2, ··· do
使用1.19:
˜
I
n−1
→
˜
I
n
end for
的误差!下面就来看看这个误差是怎么样产生的。
使用误差传递公式
10
,有
E
0
= I
0
−
˜
I
0
E
n
= I
n
−
˜
I
n
(1.21)
7
(2025.2.25) 对,这部分是算法内容,看看他讲不讲,不讲的话之后会提到的,也不做了。
8
(2025.2.27) 他讲了
9
e
−1
≈ 1 + (−1) +
(−1)
2
2!
+ · · · +
(−1)
k
k!
(1.20)
其中取不同的 k 有效数字的个数也不同。
10
也就是那个含导数用微分中值定理推出来的公式
14
数值分析笔记 目录 Tsui Dik Sang
因此
E
n
= (−1)
n
n!E
0
(1.22)
这是阶乘形式的误差!在
n
很大的时候误差的增长速度也很恐怖。
1.2.1.2 算法 B: 反向递推法
使用积分不等式找到原式的两个极值
e
−1
n + 1
= e
−1
( min
0⩽x≦1
e
x
)
ˆ
1
0
x
n
dx < I
n
< e
−1
( max
0⩽x≦1
e
x
)
ˆ
1
0
x
n
dx =
1
n + 1
(1.23)
我们从出现故障的 n=9 开始,取其近似值
11
I
9
≈
1
2
1
10
+
e
−1
10
= I
∗
9
= 0.0684 (1.24)
然后用反向的递推公式
I
∗
9
= 0.0684
I
∗
n−1
=
1
n
(1
−
I
∗
n
), n = 9, 8, ··· , 1
(1.25)
算法如下 由1.25同理可推误差
算法 2 算法 B
已知 I
∗
9
= 0.0684
for all n = 9, 8, ··· , 1 do
使用递推公式1.25:I
∗
n
→ I
∗
n−1
end for
E
0
=
(−1)
m
m!
E
∗
m
(1.26)
实际上,这个算法 B 算不上是一个通用算法,而更应该说是算法 A 的一个纠错算法,从算法 A 找出的出明显错误的地方 n=m
开始用夹逼得到一个初值,反向推回 E
∗
0
,并且在这个过程中误差越来越小。
1.2.1.3 稳定的算法
通过上面的例子,我们就可以给出定义了。
定义 1.2.1 (算法的稳定性). 一个算法如果输入数据有误差,而在计算过程中舍入误差不增长,则称算法是稳定的,否则,算
法不稳定
这里需要注意算法 B 的初值是开始逆向递推得到的初值,在上面的例子中指 I
∗
9
1.2.2 病态问题判断
定义 1.2.2 (病态问题). 若问题的解对于输入数据的微小变化非常敏感,那么称这个问题是病态问题
注意,病态问题是问题本身的问题,靠算法无能为力,像上面算法稳定性问题是算法选取不妥当导致稳定性不强的问题,
不是问题本身的问题
下面再定义一个参数来描述问题病态的程度
11
这本质上是在限定这个值,强制使得其不会犯算法 A 的错,顺带提一句算法 A 检测到该值出故障是式1.23的功劳
15
数值分析笔记 目录 Tsui Dik Sang
定义 1.2.3 (条件数). 计算函数值 f(x) 时,若 x 有微小扰动 ∆x = x −x
∗
, 定义条件数为函数值与输入数据相对误差的比值,
即
C
p
=
|
f(x)−f (x
∗
)
f(x)
|
|
∆x
x
|
≈ |
xf
′
(x)
f(x)
| (1.27)
一般认为 C
p
⩾ 10 就认为问题是病态的了其值越大,问题就越病态,
1.2.3 避免误差伤害
书上给了很多例子,但我认为其实就只有一个思想
防止“大吃小”
比如下面的几个例子
12
,都是上面的意思
•
lg x
1
− lg x
2
√
← lg
x
1
x
2
(1.28)
•
√
x + 1 −
√
x
√
←
1
√
x + 1 +
√
x
(1.29)
•
((···((10000 + δ
1
) + δ
2
) ···) + δ
n
)
√
← 10000 +
n
X
i
δ
i
(1.30)
其中 δ
i
∈ [0.00000001, 0.00000009]
1.2.4 算法设计
后面的章节就是在介绍下面的算法,基本会有重复,这里就从略了,只是稍稍提一下秦九韶算法
算法 3 秦九韶算法
输入: n + 1 个系数 a
0
, a
1
, ··· , a
n
和 x
输出: P
n
= a
0
+ a
1
x + a
2
x
2
+ ··· + a
n
x
n
1: P
1
= a
n−1
+ xa
n
2: P
2
= a
n−2
+ x(a
n−1
+ xa
n
)
.
.
.
3: P
n
= a
0
+ x(a
1
+ x(a
2
+ ··· + x(a
n−1
+ xa
n
) ···))
1.2.5 加权平均的松弛技术
定义 1.2.4 (加权平均的松弛技术). 设 x
∗
1
, x
∗
2
, ··· , x
∗
n
是 n 个近似值,ω
1
, ω
2
, ··· , ω
n
是权重系数,ω
1
+ ω
2
+ ···+ ω
n
= 1,那
么加权平均值为
x
∗
= ω
1
x
∗
1
+ ω
2
x
∗
2
+ ··· + ω
n
x
∗
n
(1.31)
其中 ω
i
的选取可以是任意的,甚至可以是负数。选取合适的 ω
i
可以使得加权平均值比原来的近似值更接近真实值。
12
“
√
←”表示右边的算法优于左边
16
第二章 插值法
2.1 插值问题与多项式插值简介
2.1.1 插值问题
先不看是不是多项式插值,我们就讨论:什么是插值
定义 2.1.1 (插值法的一系列定义). 设函数 y = f (x) 在区间 [a, b] 上有定义,且已知在点 a ⩽ x
0
< x
1
< ··· < x
n
⩽ b 上的
值分别为 y
1
, ··· , y
n
,若存在函数 P (x),使得
P (x
i
) = y
i
, i = 0, 1, ··· , n (2.1)
那么称
• P (x) 为插值函数
• 点 x
0
, x
1
, ··· , x
n
为插值点
• 包含插值节点的区间 [a, b] 为插值区间
• 求插值函数 P (x) 的过程为插值法
2.1.2 多项式插值
定义 2.1.2 (多项式插值). 若插值函数 P (x) 是一个 n 次多项式
P (x) = a
0
+ a
1
x + a
2
x
2
+ ··· + a
n
x
n
(2.2)
那么称 P (x) 为 n 次多项式插值函数
除此之外还有分段插值、三角插值等等。用自由度的思想很容易知道
1
1
严谨的证明使用线代的知识,判断插值函数是否唯一,即证明下面方程组是否有唯一解
a
0
+ a
1
x
0
+ a
2
x
2
0
+ · · · + a
n
x
n
0
= y
0
a
0
+ a
1
x
1
+ a
2
x
2
1
+ · · · + a
n
x
n
1
= y
1
· · ·
a
0
+ a
1
x
n
+ a
2
x
2
n
+ · · · + a
n
x
n
n
= y
n
(2.3)
此方程的系数矩阵为 Vandermonde 矩阵:
A =
1 x
0
x
2
0
· · · x
n
0
1 x
1
x
2
1
· · · x
n
1
· · ·
1 x
n
x
2
n
· · · x
n
n
(2.4)
17
数值分析笔记 第二章 插值法 Tsui Dik Sang
2.1.2.1 多项式插值的唯一性
定理 2.1.1 (多项式插值的唯一性). 对于有 n 个插值节点的插值问题,其有唯一的 n 次插值函数。
下面介绍的两种插值法都是多项式插值法的特例
2.2 拉格朗日插值法
2.2.1 一次插值
即先考虑多项式次数 n 为 1 的情况,此时应该有两个插值节点。我们根本不需要使用插值法,“两点确定一条直线”的结论
即可确定插值直线 (函数),
引理 2.2.1 (一次插值结果). 假设插值点为 (x
k
, y
k
), (x
k+1
, y
k+1
), 那么插值函数为
L
1
(x) =
x − x
k+1
x
k
− x
k+1
y
k
+
x − x
k
x
k+1
− x
k
y
k+1
, 两点式 (2.6)
或者表示成点斜式
L
1
(x) = y
k
+
y
k+1
− y
k
x
k+1
− x
k
(x − x
k
), 点斜式 (2.7)
我们在拉格朗日法中关注的是两点式,这个关注点与高中用向量法表示的三点共线非常相似,都是一种加权线性相加。我们也将
上面的方程写成加权线性相加的形式
L
1
(x) = y
k
l
k
(x) + y
k+1
l
k+1
(x) (2.8)
其中
l
k
(x) =
x−x
k+1
x
k
−x
k+1
l
k+1
(x) =
x−x
k
x
k+1
−x
k
(2.9)
并且其满足
l
k
(x
k
) = 1 , l
k
(x
k+1
) = 0
l
k+1
(x
k
) = 0 , l
k+1
(x
k+1
) = 1
(2.10)
定义 2.2.1 (线性插值基函数). 上面的 l
k
(x), l
k+1
(x) 称为线性插值基函数
2.2.2 Lagrange 插值公式 (多次插值)
推广到一般的情况
2
,
定义 2.2.2 (n 次插值基函数). 若 n 次多项式 l
j
(x) 在 n + 1 个节点 x
0
< x
1
< ··· < x
n
上满足
l
j
(x
i
) =
1, k = j
0, k ̸= j
j, k = 0, 1, ···n (2.11)
其行列式为
det A =
0⩽i<j⩽n
(x
j
− x
i
) ̸= 0 (2.5)
即知解 a
0
, a
1
, · · · , a
n
唯一
2
书本上其实还详细说了二次插值增进理解,这里就略过,直接给一个普遍结论了
18
数值分析笔记 第二章 插值法 Tsui Dik Sang
就称这 n + 1 个 n 次多项式为节点 x
0
, x
1
. ··· , x
n
上的 n 次插值基函数
使用零点式可以推导出
3
推论 2.2.2 (n 次插值基函数的表达式). n 次插值基函数的表达式为
l
k
(x) =
(x − x
0
) ···(x −x
k−1
)(x − x
k+1
)(x − x
n
)
(x
k
− x
0
) ···(x
k
− x
k−1
)(x
k
− x
k+1
) ···(x
k
− x
n
)
, k = 0, 1, ··· , n (2.14)
于是插值多项式 L
n
(x) 可以表示为
L
n
(x) =
n
X
k=0
y
k
l
k
(x) (2.15)
其显然是满足2.1的,即经过所有节点。
2.2.3 插值余项与误差估计
2.2.3.1 用零点式 ω 表示的插值函数
引入记号
ω
n+1
(x) = (x − x
0
)(x − x
1
) ···(x −x
n
) (2.16)
容易求得
4
ω
′
n+1
(x
k
) = (x
k
− x
0
) ···(x
k
− x
k−1
)(x
k
− x
k+1
) ···(x
k
− x
n
) (2.20)
定义 2.2.3 (余项). 若在 [a, b] 上用 L
n
(x) 近似 f (x),则可以定义截断误差,也就是所谓的余项.
R
n
(x) = f(x) − L
n
(x) (2.21)
2.2.3.2 定理阐释
定理 2.2.3 (Lagrange 插值余项). 若 f(x) ∈ C
n+1
[a, b], 且 x
0
, x
1
, ··· , x
n
为 [a, b] 上的 n+1 个节点,那么对于任意 x ∈ [a, b],
存在 ξ ∈ (a, b),使得
R
n
(x) =
f
(n+1)
(ξ)
(n + 1)!
ω
n
(x) (2.22)
3
由零点式可以写 l
k
(x) 为
l
k
(x) = A(x − x
0
) · · · (x − x
k−1
)(x − x
k+1
)(x − x
n
) (2.12)
再代入 l
k
(x
k
) = 1 可以解得待定系数 A 为
A =
1
(x
k
− x
0
) · · · (x
k
− x
k−1
)(x
k
− x
k+1
) · · · (x
k
− x
n
)
(2.13)
因此也就得证了。
4
从结果出发猜想,提出 (x − x
k
) 这一项,则
ω
n+1
(x) = (x − x
k
)k(x), k(x) = (x
k
− x
0
) · · · (x
k
− x
k−1
)(x
k
− x
k+1
) · · · (x
k
− x
n
) (2.17)
所以, 使用乘法的求导公式,有
ω
′
n+1
(
x
) =
k
(
x
) + (
x
−
x
k
)
k
′
(x) (2.18)
由 k(x
k
) = 0 得
ω
′
n+1
(x
k+1
) = k(x) = (x
k
− x
0
) · · · (x
k
− x
k−1
)(x
k
− x
k+1
) · · · (x
k
− x
n
) (2.19)
即原式证得。
19
数值分析笔记 第二章 插值法 Tsui Dik Sang
其中 ω
n
(x) 由2.20得。
并且这里一般 f
(n+1)
(ξ) < 1 定理的证明需要反复用罗尔定理。
根据这个余项定理可以验证2.1.1了
5
2.2.3.3 定理证明
显然,在节点处无误差,因此
R
n
(x
k
) = 0 , k = 0, 1, ··· , n (2.23)
于是可以用零点式来表示 R
n
(x)
R
n
(x) = K(x)(x − x
0
)(x − x
1
) ···(x −x
n
) = K(x)ω
n+1
(x) (2.24)
其中 K(x) 是待定函数。下面把 x 也看成一个定值,并做一个关于 t 的函数
φ(t) = f(t) − L
n
(t) − K(x)ω
n+1
(t) (2.25)
6
• 根据 f 的 n 阶导的连续性以及 n+1 阶导的存在性,
• 可知 φ
(n+1)
(t) 在 (a, b) 存在
• 又 φ 在 x
0
, x
1
, ··· , x
n
, x 处均为零 (有 n+2 个零点)
• 因此 φ
′
(t) 至少有 n+1 个零点
7
.
.
.
• 以此类推,φ
(n+1)
(t) 至少有一个零点 ξ ∈ (a, b), 使得
φ
(n+1)
(ξ) = f
(n+1)
(ξ) − (n + 1)!K(ξ) = 0 (2.26)
从而得到
K(x) =
f
(n+1)
(ξ)
(n + 1)!
, ξ ∈ (a, b)并且依赖于 x (2.27)
代入2.24即可得证。
实际上,这个证明和插值多项式的证明有很多类似之处:零点式,待定系数 (这里是待定函数)
2.2.3.4 余项的估计
上面精确求解余项需要知道 f
(n+1)
(ξ),这个是很难求的,因此我们需要一个估计值。
也即用 f
(n+1)
(ξ) 来估计
推论 2.2.4. 如果 max
a⩽x⩽b
|f
(n+1)
(x)| = M
n+1
, 则
|R
n
(x)| ⩽
M
n+1
(n + 1)!
|ω
n+1
(x)| (2.28)
5
即如果用 n 次多项式插值 n 次多项式,那 f
(n+1)
(ξ) = 0 , 也就意味着余项为零没有误差了。
6
前几天 (2025.2.28) 一直卡在这里,不明白为什么可以将 x 也看成定值,现在清楚了,因为 φ(t) 并不是对2.21用 t 对 x 的完全替换,K(x) 是保留了 x 的
因此 φ(t) 对所有的 t 并不是一直为零。
7
根据罗尔定理,φ
′
(t) 在 φ(t) 的两个,零点之间至少有一个零点
20
数值分析笔记 第二章 插值法 Tsui Dik Sang
2.2.3.5 余项定理的应用
这个应用是对定理2.1.1从另外一个视角的证明
利用余项定理,当 f (x) = x
k
且 k ⩽ n 时,由于 f
(n+1)
(x) = 0, 所以
R
n
∗ (x) = x
k
−
n
X
i=0
= 0 (2.29)
⇒ x
k
=
n
X
i=0
l
i
(x)x
i
(2.30)
式2.30就是说明了:用 n 个 n 次插值基函数插值 n 次多项式,得到的结果没有误差
2.3 Newton 插值法
这种插值法依据的是另外一种表示方法,如果是一次的话,即点斜式2.7。
2.3.1 插值多项式的逐次生成以及均差的引入
首先,为什么要有这种差值法呢?因为在拉格朗日插值法中,每次增加一个节点,都要重新计算一次插值函数,这样的计算
量是很大的。
8
2.3.1.1 一次差值
即只有一个点的情况,就是两点确定一条直线的情况
P
1
(x) = f(x
0
) +
f(x
1
) − f(x
0
)
x
1
− x
0
(x − x
0
) = P
0
(x) + a
1
(x − x
0
) (2.31)
其中 a
1
=
f(x
1
)−f(x
0
)
x
1
−x
0
, 从2.31可以发现,一次差值其实已经是零次差值的扩展了
2.3.1.2 二次差值
9
如果我们构造
P
2
(x) = P
1
(x) + a
2
(x − x
0
)(x − x
1
) (2.32)
就可以满足差值条件
P
2
(x
0
) = f(x
0
)
P
2
(x
1
) = f(x
1
)
P
2
(x
2
) = f(x
2
)
(2.33)
待定系数法可求得
a
2
=
P
2
(x
2
) − P
1
(x
2
)
(x
2
− x
0
)(x
2
− x
1
)
=
f(x
2
)−f(x
0
)
x
2
−x
0
−
f(x
1
)−f(x
0
)
x
1
−x
0
x
2
− x
1
(2.34)
8
当时看到这里的时候我有另外一个疑问,那 Lagrange 的优势是什么? 后面我发现,这需要回到1.2.3,Newton 插值法属于是逐个相加计算,如果后面加的
每一项都很小的话就容易出现“大吃小”的情况,而 Lagrange 法就不会
9
这部分扩展到高次有一点难理解,因此也把二次的写下来了
21
数值分析笔记 第二章 插值法 Tsui Dik Sang
2.3.1.3 均差的引入与推广到多次
如果按照上面的推法,实际就是零点式,一路推 p 下去就有
P
n
(x) = a
0
+ a
1
(x − x
0
) + a
2
(x − x
0
)(x − x
1
) + ··· + a
n
(x − x
0
)(x − x
1
) ···(x −x
n−1
) (2.35)
可是 a
n
的表示也将相当复杂,如何添加新定义来简介表明它,并且最好能让各项之间的关系也明晰呢?
定义 2.3.1 (均差). • 称 f [x
0
, x
k
] =
f(x
k
)−f(x
0
)
x
k
−x
0
为函数 f (x) 关于 x
0
, x
k
的一阶均差
• f[x
0
, x
1
.x
k
] =
f[x
0
,x
k
]−f[x
0
,x
k
]
x
k
−x
1
为函数 f (x) 关于 x
0
, x
1
, ··· , x
k
的二阶均差
• 以此类推可得 k 阶均差
f[x
0
, x
1
, ··· , x
k
] =
f[x
1
, x
2
, ··· , x
k−2
, x
k
] − f[x
0
, x
1
, ··· , x
k−1
]
x
k
− x
k−1
(2.36)
刚刚我们的系数也 a
1
, a
2
, ··· , a
n
就可以表示成 f [x
0
, x
1
], f[x
0
, x
1
, x
2
], ··· , f [x
0
, x
1
, x
2
, ···x
k
]
定理 2.3.1 (均差的线性表达形式).
f[x
0
, x
1
, ··· , x
k
] =
k
X
i=0
f(x
i
)
k
Y
j=0,j̸=i
(x
i
− x
j
)
(2.37)
上面的定理确定好 f (x
i
) 各个系数即可证明,形式上来看是分母项缺失后求和,由此有一个有用的推论
推论 2.3.2 (节点排列次序无关性). 均差与节点排列顺序无关,即
f[x
0
, x
1
, ··· , x
k
] = f[x
1
, x
0
, x
2
, ··· , x
k
] = ··· = f [x
k
, x
0
, x
1
, ··· , x
k−1
] (2.38)
将2.36换 x
k
−
1
为其他点, 比如 x
0
,就可以得到
定理 2.3.3 (均差的任一点表示). 或者可以理解为均差的递推公式
f[x
0
, x
1
, ··· , x
k
] =
f[x
1
, x
2
, ··· , x
k
] − f[x
0
, x
1
, ··· , x
k−1
]
x
k
− x
0
(2.39)
递推公式无需死记硬背,最妙的当属从形状上去理解它,看表2.1,无论是几阶均差,都可以一路回溯到金字塔的底部两个脚,对
应的就是2.39中 x
k
, x
0
,然后就可以相减了,
定理 2.3.4 (均差与导数的关系).
f[x
0
, x
1
, ··· , x
k
] =
f
(k)
(ξ)
k!
, ξ ∈ [a, b] (2.40)
上面的定义我想到的是用2.21以及多项式插值的唯一性来证明
10
,由此有一个均差形式的余项估计值
10
课本上说的罗尔定理证明不知道是不是这个意思还是另有他法
22
数值分析笔记 第二章 插值法 Tsui Dik Sang
推论 2.3.5 (均差形式的余项估计).
R
(
x) = f[x
0
, x
1
, ··· , x
n
]ω
n+1
(x) (2.41)
这个定理在已知均差表情况下的插值问题非常好用
2.3.1.4 Newton 均差插值多项式
有了前面的定义,就可以给出表示和简洁的定理了
定理 2.3.6 (Newton 均差插值多项式). 形如
P
n
(x) = a
0
+ a
1
(x − x
0
) + a
2
(x − x
0
)(x − x
1
) + ··· + a
n
(x − x
0
)(x − x
1
) ···(x −x
n−1
) (2.42)
且 a
k
= f [x
0
, x
1
, ··· , x
k
], k = 0, 1, 2, ··· , n 的称为 Newton 均差插值多项式
表 2.1: 均差表
x
i
f(x
i
) 一阶均差 二阶均差 三阶均差 四阶均差 五阶均差
x
0
f(x
0
)
f[x
0
, x
1
]
x
1
f(x
1
) f[x
0
, x
1
, x
2
]
f[x
1
, x
2
] f[x
0
, x
1
, x
2
, x
3
]
x
2
f(x
2
) f[x
1
, x
2
, x
3
] f[x
0
, x
1
, x
2
, x
3
, x
4
]
f[x
2
, x
3
] f[x
1
, x
2
, x
3
, x
4
] f[x
0
, x
1
, x
2
, x
3
, x
4
, x
5
]
x
3
f(x
3
) f[x
2
, x
3
, x
4
] f[x
1
, x
2
, x
3
, x
4
, x
5
]
f
[
x
3
, x
4
] f
[
x
2
, x
3
, x
4
, x
5
]
x
4
f(x
4
) f[x
3
, x
4
, x
5
]
f[x
4
, x
5
]
x
5
f(x
5
)
当然,有时候阶梯型的比金字塔型的更好理解,这个看个人
由于均差概念比较抽象,因此建议在之后做 Newton 或者 Hermite 插值的题的时候一般都先根据已知值画出这个表 (可
能阶数不同),方便后续计算
2.3.2 差分形式的 Newton 插值多项式
如果各节点等距,那么表示结果会进一步简化,不过在这之前,引入一个很重要的概念差分
2.3.2.1 差分算子
定义 2.3.2 (等距点). 若节点之间满足
x
k
= x
0
+ kh (2.43)
则称这一系列节点是等距点,h 被称为步长
23
数值分析笔记 第二章 插值法 Tsui Dik Sang
定义 2.3.3 (差分). 若 f
k
= f (x
k
), 则称
∆f
k
= f
k+1
− f
k
(2.44)
为 x
k
以 h 为步长的一阶差分, 类似的,称
∆
2
f
k
= ∆f
k+1
− ∆f
k
(2.45)
为 x
k
以 h 为步长的二阶差分一般的,称
∆
k
f
k
= ∆
k−1
f
k+1
− ∆
k−1
f
k
(2.46)
为 x
k
以 h 为步长的 k 阶差分
2.3.2.2 差分算子的运算
单独的差分算子 ∆ 并不能很好的进行不变换,因此再引入两个算子
定义 2.3.4 (不变算子).
I
f
k
=
f
k
(2.47)
定义 2.3.5.
Ef
k
= f
k+1
(2.48)
于是就可以用上面的两个算子来表示差分算子
推论 2.3.7 (一阶差分算子的表示).
∆f
k
= (E − I)f
k
(2.49)
进一步可以一般化
推论 2.3.8 (n 阶差分算子的表示).
∆
n
f
k
=
n
X
j=0
(−1)
j
n
j
!
f
n
+
k
−
j
(2.50)
上面的证明需要用到二项式定理
11
反之也可以用各阶差分来表示函数值 (就是一个线性变换)
推论 2.3.9 (函数值关于差分的线性表示).
f
n+k
=
n
X
j=0
n
j
!
∆
j
f
k
(2.52)
证明方法类似
12
11
∆
n
f
k
= (E − I)
n
f
k
=
n
j=0
(−1)
j
n
j
E
n−j
f
k
=
n
j=0
(−1)
j
n
j
f
n+k−j
(2.51)
12
f
n+k
= E
n
f
k
= (I + ∆)
n
f
k
==
n
j=0
n
j
∆
j
f
k
=
n
j=0
n
j
∆
j
f
k
(2.53)
24
数值分析笔记 第二章 插值法 Tsui Dik Sang
2.3.2.3 Newton 前插公式
先再给出一个均差与差分的关系
定理 2.3.10 (均差与差分).
f[x
k
, x
k+1
, ··· , x
k+m
] =
∆
m
f
k
m!h
m
(2.54)
于是代入之的插值公式就有了
定理 2.3.11 (Newton 前插公式).
P (x
0
+ th) = f
0
+ t∆f
0
+ ··· +
t(t − 1)(t − 2) ···(t − n + 1)
n!
∆
n
f
0
(2.55)
同理也容易推出其余项为
推论 2.3.12 (牛顿前插公式的余项).
R
n
(x) =
t(t − 1)(t − 2) ···(t − n + 1)
(n + 1)!
h
n+1
f
(n+1)
(ξ), ξ ∈ (x
0
, x
n
) (2.56)
2.4 Hermite 插值
之前的两种插值方法都是针对只知道函数值的情况的插值,在某一些物理场景中,可能会遇到给出了混合给函数值以及导数
的情况,此时单纯的 Newton 插值法就不够用了,这时候就需要 Hermite 插值法了。
2.4.1 函数导数的均差表示
2.4.1.1 重节点均差
引理 2.4.1 (均差的连续性). 设 f ∈ C
n
[a, b], x
0
, x
1
, ··· , x
n
为 [a, b] 上的相异节点,则 f[x
0
, x
1
, ··· , x
n
] 是关于 x
0
, x
1
, ··· , x
n
的连续函数
这个定理为下面讨论均差连续性做准备
13
定义 2.4.1 (重节点的一阶均差).
f[x
0
, x
0
] = lim
x→x
0
f[x
0
, x] = f
′
(x
0
) (2.57)
定义 2.4.2 (重节点的二阶均差).
f[x
0
, x
0
, x
1
] =
f[x
0
, x
1
] − f[x
0
, x
0
]
x
1
− x
0
(2.58)
13
这个定理的证明复杂,这里 (包括书 [2] 上) 也没有给出
25
数值分析笔记 第二章 插值法 Tsui Dik Sang
推论 2.4.2 (重节点的二阶均差与原函数的关系).
f[x
0
, x
0
, x
1
] = lim
x
1
→x
0
x
2
→x
0
f[x
0
, x
1
, x
2
] =
1
2
f
′′
(x
0
) (2.59)
这个定理的证明需要用到泰勒公式
14
。
于是,我们就可同理扩展到 n 阶
推论 2.4.3 (n 阶重节点均差).
f[x
0
, x
0
, ··· , x
0
] = lim
x
i
→x
0
f[x
0
, x
1
, x
2
, ··· , x
n
] =
1
n!
f
(n)
(x
0
) (2.60)
2.4.2 Taylor 插值多项式
2.4.2.1 多项式得到
于是对 Newton 均差插值多项式2.42进行取极限操作
推论 2.4.4 (Taylor(插值) 多项式).
P
n
(x) = lim
x
i
→x
0
P
n
(x) = f(x
0
) + f
′
(x
0
)(x − x
0
) + ··· +
f
(n)
(x
0
)
n!
(x − x
0
)
n
i = 1 , 2, 3, ··· , n (2.61)
容易看到,式2.61满足含导数值的插值条件
P
(k)
n
(x
0
) = f
(k)
(x
0
) (2.62)
因此可以用这个多项式来解决本节开头提到的问题。
2.4.2.2 余项误差
对于完全只给定一个点的 n 重导数值的插值问题——这本质上就是一个泰勒展开!那么余项误差的公式就是泰勒展开余项
公式,
推论 2.4.5 (Taylor 插值的余项).
R
n
(x) =
f
(n+1)
(ξ)
(n + 1)!
(x − x
0
)
n+1
, ξ ∈ (a, b) (2.63)
2.4.2.3 一般规律总结
上面的是只给导数的理想情况,对于一般,我们其实不是需要几个导数或者几个值,而是以条件来衡量
定理 2.4.6 (Hermite 插值适用条件). 只要给出独立的m + 1 个插值条件 (混合含函数值与导数值),就可以得到次数不超过 m
次的 Hermite 插值多项式。如果不独立,则需要具体分析条件是否矛盾
对于这个定理,一定要搞清楚给的条件是否独立,在期末复习的时候我看到了一道很有意思的题,
14
将 f (x) 在 x
0
处的泰勒展开代入式2.58,(后面的项用高阶小项符号表示,之后可约),然后计算极限即可。
26
数值分析笔记 第二章 插值法 Tsui Dik Sang
1. 求一个 3 次多项式使之满足
H
3
(a) = f(a), H
′
3
(a) = f
′
(a), H
3
(b) = f(b), H
′
3
(b) = f
′
(b), (2.64)
2. 讨论是否存在一个 4 次多项式 (也就是说其四次项不为零),使得
H
3
(a) = f(a), H
′
3
(a) = f
′
(a), H
3
(b) = f(b), H
′
3
(b) = f
′
(b), H
′
a + b
2
= f
′
a + b
2
(2.65)
第一道是常规题,用下面说的2.4.3.2可以直观清晰的求解
15
然后第二题在我没有意识到独立之前,我一眼认为是存在的,实
际上,答案是不一定,虽然是有五个条件,但是第五个条件不独立!何以见得?因为可以用前面四个条件推出第五个条件,如果
矛盾的话就不存在的,不矛盾的话就可能存在,题目的证明方法也正是如此
16
推论 2.4.7 (Hermite 插值次数). 能确定的次数应该比条件数小 1。
2.4.3 Hermite 插值插值:利用重节点均差
由于是混合,情况将非常多样,
需要具体问题具体分析,给通式只会让人觉得晕,所以我还是以一道例题切入,忘读者能悟出思路
求一个次数不高于 4 的多项式 P (x), 使得
P
4
(1) = 2, P
4
(2) = 4, P
4
(3) = 12, P
′
4
(1) = 1, P
′
4
(3) = −1 (2.73)
15
当然,计算量肯定会很难受的了,但是至少能保证思路是清晰的
16
首先由于两题给的条件前四个都有,所以 H(x), H
3
(x) 均满足前四个条件,写成因式形式如下
H(x) − H
3
(x) = A(x − a)
2
(x − b)
2
H
′
(x) − H
′
3
(x) = aA(x − a)(x − b)(2λ − a − b)
(2.66)
然后我们代入第五个条件的点得到
H
′
a + b
2
− H
′
3
a + b
2
= 0 (2.67)
所以
H
′
3
1 + b
2
= H
′
1 + b
2
(2.68)
这意味着如果 H(x) 满足的这个第五个条件 H
3
(x) 也满足,而 H
3
(x) 我们已经可以通过前四个条件确定下来了,
H
′
3
(x) = f
′
(a)x − 2f [a, a, b](x − a) + f [a, a, b, b](x − a)(3x − 2a − b) (2.69)
⇒ H
′
3
a + b
2
= f
′
(a) + f [a, a, b](b − a) − f [a, a, b, b]
(b − a)
4
4
(2.70)
那么如果
•
f
′
(a) + f [a, a, b](b − a) − f [a, a, b, b]
(b − a)
4
4
= f
′
a + b
2
(2.71)
即第五个条件不矛盾,那么 H(x) 存在,其与 H
3
(x) 满足关系
H(x) = H
3
(x) + A(x − a)
2
(x − b)
2
, A ∈ R (2.72)
• 否则,第五个条件矛盾,不存在这样的四次多项式 H(x)
也就是说 H
(
x) 就算不给第五个条件也已经确定了,
27
数值分析笔记 第二章 插值法 Tsui Dik Sang
2.4.3.1 待定系数法
由2.4.6可以知道确定的多项式是四次,那么就设
P
4
(x) = a
0
+ a
1
x + a
2
x
2
+ a
3
x
3
+ a
4
x
4
(2.74)
求导数,然后根据上面的五个条件列出方程组,求这五个系数就可以了啊,优点是思路不容易错,缺点是计算可能容易错
2.4.3.2 重节点均差法
我们设
x
0
= 1
x
1
= 2
x
2
= 3
,于是多项式就可以表示成
17
P (x) = f(x
0
)+f[ x
0
, x
0
](x−x
0
)+f[x
0
, x
0
, x
1
](x−x
0
)
2
+f[x
0
, x
0
, x
1
, x
2
](x−x
0
)
2
(x−x
1
)+f[x
0
, x
0
, x
1
, x
2
, x
2
](x−x
0
)
2
(x−x
1
)(x−x
2
)
(2.75)
然后列出形如2.1的表格,根绝表格数据填入多项式,就可得出结果了。
这个方法书上没有提到,但是我是最好最容易理解思路最清晰的方法!
2.4.4 两个典型的 Hermite 插值
但是我必须很遗憾的告诉你,这是书上的方法,后来在我做题的过程中被证明及其难用,但是当时做笔记的时候以为是唯
一方法记录了下来,所以这一部分实际上可以跳过,掌握上面的 Hermite 插值的通式就可以了。
2.4.4.1 P (x
i
) = f(x
i
)(i = 0, 1, 2) & P
′
(x
1
) = f
′
(x
1
)
由定理2.4.6知道对于这个条件需要 3 次插值多项式。
• 先由给定的函数值确定出一个待定系数的形式
P (x) = f(x
0
) + f[x
0
, x
1
](x − x
0
) + f[x
0
, x
1
, x
2
](x − x
0
)(x − x
1
) + A(x − x
0
)(x − x
1
)(x − x
2
) (2.76)
前面两项是为插值点个数而定的二次项,而最后一项是为了被导数值限定而增添的待定系数项。通过导数值条件代入可以
解得
A =
f
′
(x
1
) − f[x
0
, x
1
] − (x
1
− x
0
)f[x
0
, x
1
, x
2
]
(x
1
− x
0
)(x
1
− x
2
)
(2.77)
代入就得到了插值函数
• 然后还需要求其余项公式:
定理 2.4.8 (Hermite 插值的余项).
R(x) =
1
4!
f
(4)
(ξ)(x − x
0
)
2
(x − x
1
)
2
(2.78)
证明方法与2.2.3类似,即构造 φ(t)
17
构造次序可以有所不同,直观理解的规则就是:均差中出现的新数将在下一项的因式中出现
28
数值分析笔记 第二章 插值法 Tsui Dik Sang
2.4.4.2 两点三次多项式插值
即满足差值条件
H
3
(x
k
) = y
k
, H
3
(x
k+1
) = y
k+1
H
′
3
(x
k
) = m
k
, H
′
3
(x
k+1
) = m
k+1
(2.79)
下面用基函数的方法来讨论插值,另
H
3
(x) = α
k
(x)y
k
+ α
k+1
(x)y
k+1
+ β
k
(x)m
k
+ β
k+1
(x)m
k+1
(2.80)
其中 α
k
, α
k+1
函数分别就是关于函数值的基函数,而 β
k
, β
k+1
函数分别就是关于导数值的基函数
18
。然后通过三次函数的待定
系数法
19
, 就可以算出四个基函数的形式
α
k
(x) =
1 + 2 ·
x − x
k
x
k+1
− x
k
x − x
k+1
x
k
− x
k+1
2
α
k+1
(x) =
1 + 2 ·
x − x
k+1
x
k
− x
k+1
x − x
k
x
k+1
− x
k
2
(2.84)
β
k
(x) = (x − x
k
)
x − x
k+1
x
k
− x
k+1
2
β
k+1
(x) = (x − x
k+1
)
x − x
k
x
k+1
− x
k
2
(2.85)
代入就可以得到插值多项式 (真的是太复杂了,不写了,之后具体情况具体分析)
关于余项,也是类似的证明方法得
20
R
3
(x) =
f
(4)
(ξ)
4!
(x − x
k
)
2
(x − x
k+1
)
2
.ξ ∈ (x
k
, x
k+1
) (2.86)
2.5 分段低次插值
不应该是越高次拟合效果越好吗?为什么要用低次的呢?有两个原因:
1. 次数越高,计算复杂度也越高,如果插值点非常多的话很麻烦
2. 高次插值在本质上存在无法克服的弊病!(下面会讲)
2.5.1 高次插值的病态性质
总之,就只需要知道
18
即满足
α
k
(x
k
) = 1, α
k
(x
k+1
) = 0, α
′
k
(x
k
) = α
′
k
(x
k+1
) = 0
α
k+1
(x
k
) = 0, α
k+1
(x
k+1
) = 1, α
′
k+1
(x
k
) = α
′
k+1
(x
k+1
) = 0
β
k
(x
k
) = β
k
(x
k+1
) = 0, β
′
k
(x
k
) = 1, β
′
k
(x
k+1
) = 0
β
k+1
(x
k
) = β
k+1
(x
k+1
) = 0, β
′
k+1
(x
k
) = 0, β
′
k+1
(x
k+1
) = 1
(2.81)
19
还是相当繁琐的,不过还是有一定的技巧,比如课本中就设置了 α
k
函数的待定系数为
α
k
(x) = (ax + b)
x − x
k+1
x
k
− x
k+1
2
(2.82)
β
k
(x) = a(x − x
k
)
x − x
k+1
x
k
− x
k
+1
(2.83)
20
之后的余项有敏感数字的都会标红,既然是开卷,就要在遇到相似的时候快速知道误差什么的是哪一个数字从而快速查找
29
数值分析笔记 第二章 插值法 Tsui Dik Sang
定理 2.5.1 (高次插值的病态性质). 对于任意的的插值节点,当 n → ∞ 时,L
n
(x) 不一定收敛于 f (x),
对于这个定理不需要知道证明,知道 Runge 给出的例子 f(x) =
1
1+x
2
在 [-5,5] 上的插值异常就行了。因此我们需要分段线性插
值
2.5.2 分段线性插值
没什么好说的,就是用一段段折线去拟合。
定义 2.5.1 (分段线性插值). 设已知节点 a < x
0
< x
1
< ··· < x
n
= b 上的函数值分别为 f
0
, f
1
, ··· , f
n
, 记 h
k
= x
k+1
−x
k
, h =
max h
k
, 若一折线函数满足
• I
k
(x ∈ C[a, b])
• I
k
(x
k
) = f
k
, k = 0, 1, ··· , n
• I
k
(x) 在 [x
k
, x
k+1
] 上是线性函数
则称 I
k
(x) 为 f (x) 在 [x
k
, x
k+1
] 上的分段线性插值函数
根据余项公式的一次形式有
推论 2.5.2 (分段线性插值的余项).
max
x
k
⩽x⩽x
k+1
|f(x) − I
h
(x)| ⩽
M
2
2
max
x
k
⩽x⩽x
k+1
|(x − x
k
)(x − x
k+1
)| =
M
2
8
h
2
(2.87)
由此也就证明了这个方法收敛性
推论 2.5.3 (分段线性插值的收敛性).
lim
h→0
I
h
(x) = f(x) (2.88)
2.5.3 分段三次插值
其实引入三次的原因不是因为更逼近,而是为了解决满足 Hermite 插值条件的问题——即已知导数!
定义 2.5.2 (分段三次 Hermite 插值). 若在节点 x
k
, (k = 0, 1, 2, ··· , n) 上除了已知函数值 f
k
还给出了导数值 f
′
k
, (k =
0, 1, 2, ··· , n),那么构造的插值函数 I
h
(x) 满足
• I
h
(x) ∈ C[a, b]
• I
h
(x
k
) = f
k
, , I
′
h
(h) = f
′
k
, k = 0, 1, ··· , n
• I
h
(x) 在 [x
k
, x
k+1
] 上是三次函数
同样根据余项公式的三次形式有
30
数值分析笔记 第二章 插值法 Tsui Dik Sang
推论 2.5.4 (分段三次插值的余项).
max
a⩽x⩽b
= |f (x) −I
h
(x)| ⩽
M
4
384
h
4
(2.89)
2.6 三次样条插值
也就是再增加了二阶导数限制 (或者在一阶导数连续), 但是没想到他的方法会这么的繁琐抽象。
定义 2.6.1 (三次样条插值). 若函数 S(x) ∈ C
2
[a, b],
• 在每一个小区间上是三次多项式
• 满足插值条件 S(x
i
) = y
i
• 二阶导连续
则称 S(x) 为三次样条插值函数
2.6.1 条件的确定
2.6.1.1 罗列已知条件
目前已有的条件有
•
S(x
j
− 0) = S(x
j
+ 0)
S
′
(x
j
− 0) = S
′
(x
j
+ 0)
S
′′
(x
j
− 0) = S
′′
(x
j
+ 0)
j = 1, 2, ··· , n − 1 (2.90)
有 3n-3 个条件
• n+1 个插值点又有 n+1 个条件
总共有 4n-2 个,
2.6.1.2 对比目标所需条件数
而题目的 S(x) 在 n 个区间里都是三次函数,也就是每个区间要四个系数,总共需要 4n 个条件。
2.6.1.3 尝试构造缺失的条件
因此现在需要通过边界值的方式再增加两个条件!
• 增加一阶导的边界
S
′
(x
0
) = f
′
0
S
′
(x
n
) = f
′
n
(2.91)
• 增加二阶导边界
S
′′
(x
0
) = f
′′
0
S
′′
(x
n
) = f
′′
n
(2.92)
当 f
′′
0
= f
′′
n
= 0 的时候称为自然边界条件。
31
数值分析笔记 第二章 插值法 Tsui Dik Sang
• 给周期边界条件
S(x
0
− 0) = S(x
n
− 0)
S
′
(x
0
+ 0) = S
′
(x
n
− 0)
S
′′
(x
0
+ 0) = S
′′
(x
n
− 0)
(2.93)
2.6.2 样条插值函数的建立
三次函数不好待定系数,就用其二阶导 (也就是一次函数了) 来待定, 设
S
′′
(x) = M
j
x
j+1
− x
h
j
+ M
j−1
x − x
j
h
j
(2.94)
然后经历一系列复杂的待定系数方程
21
最后得到方程组
µ
j
M
j−1
+ +aM
j
+ λ
j
M
j+1
= d
j
µ
j
=
h
j−1
h
j−1
+h
j
λ
j
=
h
j
h
j−1
+h
j
d
j
= 6 ·
f[x
j
,x
j+1
]−f[x
j−1
,x
j
]
h
j−1
+h
j
= 6f [x
j−1
x
j
, x
j+1
]
j = 1, 2, ··· , n − 1 (2.95)
2.6.2.1 第一类边界条件:已知边界一阶导
不再说证明了,给结论即可:
如果令
λ
0
= 1
d
0
=
6
h
0
(f[x
0
, x
1
] − f
′
0
)
µ
n
= 1
d
n
=
6
h
n−1
(f
′
n
− f[x
n−1
, x
n
])
那么可以将2.95化成矩阵形式
22
1 λ
0
µ
1
2 λ
1
.
.
.
.
.
.
.
.
.
µ
n−1
2 λ
n−1
µ
n
2
M
0
M
1
.
.
.
M
n−1
M
n
=
d
0
d
1
.
.
.
d
n−1
d
n
(2.96)
然后一切交给计算机!
2.6.2.2 第二类边界条件:已知边界二阶导
首先在端点处直接得
M
0
= f
′′
0
, M
n
= f
′′
n
(2.97)
然后令
λ
0
= µ
0
= 0
d
0
= 2f
′′
0
d
n
= 2f
′′
n
然后也可以写出2.96的形式。
21
属实是很复杂,感兴趣的自己看书就好,这里不放了
22
没有写的元素是 0
32
数值分析笔记 第二章 插值法 Tsui Dik Sang
2.6.2.3 第三种边界条件:周期
不知道具体值,但是有关于边界的方程
M
0
= M
n
λ
n
M
n+1
+ µ
n
M
n−1
+ 2M
n
= d
n
λ
n
=
h
0
h
n−1
+h
0
µ
n
= 1 − λ
n
=
h
n−1
h
n−1
+h
0
d
n
= 6 ·
f[x
0
,x
1
]−f[x
n−1
,x
n
]
h
0
+h
n−1
(2.98)
写成矩阵形式:
2 λ
1
µ
1
µ
2
2 λ
2
.
.
.
.
.
.
.
.
.
µ
n−1
2 λ
n−1
λ
n
µ
n
2
M
0
M
1
.
.
.
M
n−1
M
n
=
d
0
d
1
.
.
.
d
n−1
d
n
(2.99)
实际上,记住后面这些矩阵方程以及各个字母表示什么即可了,对于开卷考试的工科同学来说足以解决实际问题了!
33
数值分析笔记 第二章 插值法 Tsui Dik Sang
34
第三章 线性方程组解法
3.1 引言:一些基本的矩阵概念
先补充一些这里用到的矩阵知识而线性代数没学或者是学过忘了的知识。读者亦可先不看,在该章之后的学习中碰到了再回
来了解。
3.1.1 非奇异矩阵
推论 3.1.1 (矩阵的行列式).
det(AB) = det(A) det(B)
det(cA) = c
n
det(A)
(3.1)
然后非奇异矩阵定义为
定义 3.1.1 (非奇异矩阵). 若矩阵 A 的行列式 det(A) ̸= 0,则称 A 为非奇异矩阵
3.1.2 特征值与谱半径
特征值的定义参见矩阵分析的笔记,而谱与半径定义为
定义 3.1.2 (谱与谱半径). 设 A 是一个 n ×n 的矩阵,λ
1
, λ
2
, ··· , λ
n
是 A 的特征值,则称 λ
1
, λ
2
, ··· , λ
n
的集合为 A 的谱,
而 A 的谱半径定义为
ρ(A) = max{|λ
1
|, |λ
2
|, ··· , |λ
n
|} (3.2)
然后需要记住三个重要结论
推论 3.1.2 (转置矩阵的特征值). λ 是 A 的特征值,那么 A 的转置矩阵 A
T
的特征值也是 λ
推论 3.1.3 (逆矩阵的特征值). λ 是 A 的特征值,那么 A 的逆矩阵 A
−1
的特征值是
1
λ
推论 3.1.4 (相似矩阵的特征值). B = S
−1
AS 是 A 的相似矩阵,那么 A 和 B 的特征值相同
3.1.3 Jordan 标准型
参见矩阵分析部分的笔记,这里不再赘述。
35
数值分析笔记 第三章 线性方程组解法 Tsui Dik Sang
3.1.4 范数
首先,什么是范数
1
,是一种“大小”的度量,是对向量“模”这个概念的引申。
3.1.4.1 向量的范数
所以向量的 p-范数定义为
定义 3.1.3 (向量的范数). 设 x = (x
1
, x
2
, ··· , x
n
) 是一个 n 维向量,x 的 p-范数定义为
∥x∥
p
=
n
X
k=1
x
p
k
!
1
p
(3.3)
其 2-范数就是我们熟悉的“模”,其他两个特殊的还有
推论 3.1.5 (1-范数与无穷范数).
∥x∥
1
=
n
X
k=1
|x
k
|,
∥x∥
∞
= max
1⩽k⩽n
|x
k
|
(3.4)
对于无穷范数的证明,取极限即可
2
,
到底哪一种范数才能表示真正的“大小”呢?对此有范数的等价性
定理 3.1.6 (向量范数的等价性). 对于一个 n 维向量 x,存在一个常数 c
1
, c
2
使得
c
1
∥x∥
p
⩽ ∥x∥
q
⩽ c
2
∥x∥
p
(3.6)
注意,这个定理不能扩展到无穷维,即 p, q 不能取到无穷
3.1.5 矩阵范数
定义 3.1.4 (矩阵的范数). 设 A 是一个 m × n 的矩阵,A 的范数定义为
∥A∥ = max
x̸=0
∥Ax∥
∥x∥
(3.7)
对于这种范数,仍然满足相容条件,即
定理
3.1.7 (
矩阵范数的相容条件
).
∥AB∥ ⩽ ∥A∥∥B∥ (3.8)
同样矩阵的范数有多种,最常见的还是那三种
1
没办法,学数值分析学到范数的时候矩阵分析还没有学到范数
2
比较小的数经过无穷次方后与最大那个数的无穷次方相比根本不在一个数量级,于是即可证明
lim
p→∞
n
k=1
|x
k
|
p
1
p
= lim
p→∞
max
1⩽k⩽n
|x
k
|
p
1
p
= lim
p→∞
max
1⩽k⩽n
|x
k
| = (3.5)
36
数值分析笔记 第三章 线性方程组解法 Tsui Dik Sang
推论 3.1.8 (矩阵的无穷范数 (行范数)).
∥A∥
∞
= max
1⩽i⩽m
n
X
j=1
|a
ij
| (3.9)
推论 3.1.9 (矩阵的 1-范数 (列范数)).
∥A∥
1
= max
1⩽j⩽n
m
X
i=1
|a
ij
| (3.10)
推论 3.1.10 (矩阵的 2-范数).
∥A∥
2
=
q
ρ(A
T
A) (3.11)
其中 ρ(A) 是矩阵的谱半径
3
,
3.1.6 向量、矩阵序列的极限
3.1.6.1 基础定义和结论
定义 3.1.6 (向量、矩阵序列的极限). 设 {x
k
} 是一个向量序列,{A
k
} 是一个矩阵序列,如果存在向量 x 和矩阵 A 使得
lim
k→∞
∥x
k
− x∥ = 0, lim
k→∞
∥A
k
− A∥ = 0 (3.12)
则称 x 为 {x
k
} 的极限,A 为 {A
k
} 的极限,记作
lim
k→∞
x
k
= x, lim
k→∞
A
k
= A (3.13)
定理 3.1.11 (向量、矩阵序列极限存在的充要条件). 设 {x
k
} 是一个向量序列,{A
k
} 是一个矩阵序列,则 {x
k
} 收敛的充分
必要条件是 {x
k
} 的每一个分量序列都收敛,即对序列的第 k 个 A
(k)
的每一个元素 a
(k)
ij
lim
k→∞
a
(k)
ij
= a
ij
, ∀i, j (3.14)
推论 3.1.12 (极限为零矩阵的线性方程组). lim A
k
= 0 的充要条件是 lim A
k
x = 0
3.1.6.2 极限与谱半径
同时,极限与谱半径之间也有一些结论
3
定义 3.1.5 (矩阵的谱半径). 设 A 是一个 n × n 的矩阵,λ
1
, λ
2
, · · · , λ
n
是 A 的特征值,
37
数值分析笔记 第三章 线性方程组解法 Tsui Dik Sang
推论 3.1.13 (三个等价的命题). 以下三个命题是等价的:
• lim
k→∞
A
k
= 0
• ρ(A < 1)
• 至少有一种从属的矩阵范数 ∥ · ∥ 使 ∥B∥
ε
< 1
推论 3.1.14 (极限等式).
lim
k→∞
∥B
k
∥
1
k
= ρ(B) (3.15)
3.1.7 对角占优矩阵等矩阵
3.1.7.1 定义
定义 3.1.7 (严格对角占优矩阵). 若矩阵 A 是一个 n × n 的矩阵,且满足
|a
ii
| >
n
X
j=1
j̸=i
|a
ij
|, i = 1, 2, ··· , n (3.16)
则称 A 是一个严格对角占优矩阵
定义 3.1.8 (弱对角占优矩阵). 若矩阵 A 是一个 n × n 的矩阵,且满足
|a
ii
| ⩾
n
X
j=1
j̸=i
|a
ij
|, i = 1, 2, ··· , n (3.17)
定义 3.1.9 (可约与不可约矩阵). 若矩阵 A 可以分解为两个矩阵 B 和 C 的直和,即
A =
A
11
A
12
0 A
22
!
(3.18)
则称 A 是可约矩阵,否则称为不可约矩阵
定义 3.1.10 (正定矩阵). 若矩阵 A 是一个 n × n 的矩阵,且满足
x
T
Ax > 0, ∀x ̸= 0 (3.19)
则称 A 是一个正定矩阵
3.1.7.2 定理
38
数值分析笔记 第三章 线性方程组解法 Tsui Dik Sang
定理 3.1.15 (对角占优定理). 若 A 是一个严格对角占优矩阵或者不可约弱对角占优矩阵,则 A 是非奇异矩阵
定理 3.1.16 (正定矩阵的特征值). 若 A 是一个正定矩阵,则 A 的特征值都是正实数
3.2 直接解法
3.2.1 高斯消元法
毫无技术含量,就是通过转换增广矩阵为上三角矩阵,然后得解,而且这是非常死板的只看首元素,所以对于 A =
0 1
1 0
!
这样的矩阵这个方法立刻就会失效
3.2.1.1 列主元高斯消元法
在每一列消元成下三角矩阵的时候多了一步:将通过行变换,让 a
kk
是最大的,从而减少由于小主元造成的在除法上的误差
这里的建议是不要记公式,公式太抽象了,记住思想就好了 这一步衍生出了一种分解方法
4
:PLU 分解
定义 3.2.1 (PLU 分解). 对于一个矩阵 A,可以分解为
A = P LU (3.20)
其中 P 为排列矩阵,L 和 U 是形如 LU 分解3.2.1的矩阵,上面的分解也可以写成
P A = LU (3.21)
的形式,注意,两个式子的 P 是不一样的
这种方法实际上是将一些可能主元为零情况下的 A 在广义上也进行了 LU 分解,下面是 ai 给我的一个 LU 分解的算法
5
算法 4 PLU 分解
输入: A
输出: 分解形式 A = P LU(或者是另外的形式)
1: 初始化
A → U
I → P , I → L
2: for all 每一列的消元 || 还没消到最后一列 do
3: 通过 (P
.
.
. U) 记录行变换,使得 U 的主元是最大
4: 通过 (L
.
.
. U) 记录分解过程,
5: end for
6: 输出 P, L, U
4
书上并没有很清晰的给出定义
5
这个算法的原理也很容易理解,两个联立矩阵的过程代表是求出 XA = U 的过程,根据变换后的结果 L, P ,即 B = L
−1
P
−1
,于是就很显然的推出了
PLU 分解的式子
39
数值分析笔记 第三章 线性方程组解法 Tsui Dik Sang
3.2.2 三角分解法
3.2.2.1 LU 分解
定理 3.2.1 (LU 分解). 对于一个非奇异矩阵 A,可以分解为一个下三角矩阵 L 和一个上三角矩阵 U 的乘积
A = LU =
1 0 ··· 0
l
21
1 ··· 0
.
.
.
.
.
.
.
.
.
.
.
.
l
n1
l
n2
··· 1
u
11
u
12
··· u
1n
0 u
22
··· u
2n
.
.
.
.
.
.
.
.
.
.
.
.
0 0 ··· u
nn
(3.22)
关于 l 元素与 u 元素与 a 元素的关系有公式,但是还是那句话,不要记公式, 具体问题具体分析,待定系数法比套公式简单 分解
之后按下面步骤求解
• Ux = L
−1
b
• x = U
−1
L
−1
b
实际上,对于计算机才会使用逆矩阵做,如果是手算的话对于上下三角矩阵式很好解的,分两步解开即可。
3.2.2.2 平方根法
是 LU 分解的扩展版,在 A 是对称阵的情况下,使用下面的分解:
定理 3.2.2 (对称阵的三角分解定理). 对于一个 n 阶对称阵 A,可以唯一分解为
A = LDL
T
=
l
11
0 ··· 0
l
21
l
22
··· 0
.
.
.
.
.
.
.
.
.
.
.
.
l
n1
l
n2
··· l
nn
d
1
0 ··· 0
0 d
2
··· 0
.
.
.
.
.
.
.
.
.
.
.
.
0 0 ··· d
n
l
11
l
21
··· l
n1
0 l
22
··· l
n2
.
.
.
.
.
.
.
.
.
.
.
.
0 0 ··· l
nn
(3.23)
分解完之后的步骤类似。但是对角阵的出现使得求逆更简单了, 这是改进的平方法,如果不改进对 A 的分解就是
A = LL
T
(3.24)
3.2.2.3 追赶法
用于处理 A 是对角占优的三对角矩阵的情况,可以将其做一下分解
定理 3.2.3 (对角占优三对角矩阵的分解).
A =
b
1
c
1
a
2
b
2
c
2
.
.
.
.
.
.
.
.
.
a
n−1
b
n−1
c
n−1
a
n
b
n
=
α
1
r
2
α
2
.
.
.
.
.
.
r
n
α
n−1
1 β
1
1 β
2
.
.
.
.
.
.
1
(3.25)
当然,在实际使用的时候还是使用待定系数法比记公式要简单。
由于这个算法的计算流程是
β
1
→ β
2
→ ··· → β
n
y
1
→ y
2
→ ··· → y
n
x
n
→ x
n−1
→ ··· → x
1
, 上面两行称为追,下面两行称为赶,因此称为追赶法。
40
数值分析笔记 第三章 线性方程组解法 Tsui Dik Sang
3.2.3 误差分析
这是对问题本身是否“病态”以及病态程度的一个度量,即系数矩阵 A 或 b 的微小变化对解的影响程度。
6
3.2.3.1 b 的变换对解的影响
定理 3.2.4 (b 的变换对解的影响). 设 A 是一个非奇异矩阵,b 是一个向量,x 是 Ax = b 的解,∆b 是 b 的微小变化,∆ x
是 x 的微小变化,那么有
∥∆x∥
∥x∥
⩽
∥A
−1
∥∥A∥∥∆b∥
∥b∥
(3.26)
这个结论经过简单的变换可以证出来
7
3.2.3.2 A 的变换对解的影响
定理 3.2.5 (A 的变换对解的影响). 设 A 是一个非奇异矩阵,b 是一个向量,x 是 Ax = b 的解,∆A 是 A 的微小变化,
∆x 是 x 的微小变化,那么有
∥∆x∥
∥x∥
⩽
∥A
−1
∥∥A∥
∥∆A∥
∥A∥
1 − ∥A
−1
∥∥A∥
∥∆A∥
∥A∥
(3.29)
证明是相似的,从略了。
3.2.4 条件数
比较定理3.2.4和3.2.5, 可以发现 ∥A
−1
∥∥A∥ 的大小与病态程度有正的关系,因此定义
定义 3.2.2 (条件数). 设 A 是一个非奇异矩阵,κ(A) 是 A 的条件数,定义为
cond(A)
v
= ∥A
−1
∥
v
∥A∥
v
(3.30)
2-范数与无穷范数产生的条件数是比较常实用的,2-范数形成的条件数又称为谱条件数,
推论 3.2.6 (谱条件数).
cond(A)
2
= ∥A
−1
∥
2
∥A∥
2
=
s
λ
max
(A
T
A)
λ
min
(AA
T
)
(3.31)
当 A 为对称矩阵的时候,
6
这里需要反复用到矩阵范数的只是,参见3.1.5
7
由
A(x + ∆x) = b + ∆b ⇒ ∆x = A
−1
∆b ⇒ ∥∆x∥ ⩽ ∥A
−1
∥∥∆b∥ (3.27)
而对线性方程组使用相容性定理3.1.7得到
∥b∥ ⩽ ∥A∥∥x∥ ⇒
1
∥x∥
⩽
∥A∥
∥b∥
b ̸= 0 (3.28)
上面推出来的两个不等式相乘就可以得到结论
41
数值分析笔记 第三章 线性方程组解法 Tsui Dik Sang
推论 3.2.7 (对称矩阵的条件数).
cond(A)
2
=
|λ
max
(A)|
|λ
min
(A)|
(3.32)
之后也没有什么好说的了。
3.3 迭代法
3.3.1 迭代法的思想
3.3.1.1 算法
首先要说的只能算是一个猜测 (算法3.3.1.1),至于 A 怎么样去拆就是各种迭代法的优劣竞争了。
算法 5 矩阵迭代法
输入: 线性方程组 Ax = b;
初值 x
(0)
;
迭代次数 k
max
或者精度 ε
输出: 近似解 x
(k)
1: 将 A 拆开:A = M − N , 方程变为 Mx = b + Nx
2: 令
M
−1
b = f
M
−1
N = B
3: while k ⩽ k
max
(或者对于限定精度的:∥x
(k+1)
− x
(k)
∥ < ε) do
4: (迭代):x
(k+1)
= Bx
(k)
+ f
5: end while
6: return x
(k)
3.3.1.2 收敛性
也可以理解为可行性, 根据定理3.1.13可以证明下面由谱半径判断的充要条件
定理 3.3.1 (谱半径判断收敛的充要条件).
ρ(B) < 1 ↔ 迭代法收敛 (3.33)
3.3.2 Jacobi 迭代法
3.3.2.1 基本定义
定理
3.3.2 (
Jacobi
迭代法
). 对于线性方程组
Ax
=
b
,可以将
A
分解为
A = D − L − U (3.34)
其中 D 是 A 的对角元素,L 是 A 的下三角元素的负数,U 是 A 的上三角元素的负数,那么可以得到迭代格式
x
(k+1)
= D
−1
(b + (L + U)x
(k)
) (3.35)
42
数值分析笔记 第三章 线性方程组解法 Tsui Dik Sang
或者说
B = D
−1
(L + U), f = D
−1
b (3.36)
通过待定系数法可以写出 Jacobi 单个元素的递推形式
8
推论 3.3.3 (Jacobi 迭代法的递推方程).
x
(0)
= (x
(0)
1
, x
(0)
2
, ··· , x
(0)
n
)
x
(k+1)
i
=
1
a
ii
b
i
−
n
X
j=1,j̸=i
a
ij
x
(k)
j
, i = 1, 2, ··· , n
(3.37)
3.3.2.2 收敛性
与 Gauss Seidel 一起讲,这里给出特征多项式
det(λI − J) = det(λI − D
−1
(L + U))
= det(D
−1
) det(λD − L − U )
=
1
det D
det C
j
(3.38)
于是,
推论 3.3.4 (Jacobi 迭代法的特征多项式). 特征多项式等价于 C
j
的行列式为零
C
j
= λD − L − U =
λa
11
a
12
··· a
1n
a
21
λa
22
··· a
2n
.
.
.
.
.
.
.
.
.
.
.
.
a
n1
a
n2
··· λa
nn
(3.39)
由于 det D ̸= 0 则算 ρ(B) 时候可以直接由 det C
j
= 0 来算
3.3.3 Gauss Seidel 迭代法
3.3.3.1 基本定义
定理 3.3.5 (Gauss Seidel 迭代法). 对于线性方程组 Ax = b,可以将 A 分解为
A = D − L − U (3.40)
其中 D 是 A 的对角元素,L 是 A 的下三角元素的负数,U 是 A 的上三角元素的负数,那么可以得到迭代格式
x
(k+1)
= (D − L)
−1
(b + Ux
(k)
) (3.41)
或者说
B = (D − L)
−1
U, f = ( D − L)
−1
b (3.42)
同样的也可以写出 Gauss Seidel 单个元素的递推形式
8
不过还是那句话,真正编程的时候公式兴许有用,但是如果是考试中矩阵阶数并不高的话直接用待定系数法是最直观最不容易出错的。
43
数值分析笔记 第三章 线性方程组解法 Tsui Dik Sang
推论 3.3.6 (Gauss Seidel 迭代法的递推方程).
x
(0)
= (x
(0)
1
, x
(0)
2
, ··· , x
(0)
n
)
x
(k+1)
i
=
1
a
ii
b
i
−
i−1
X
j=1
a
ij
x
(k+1)
j
−
n
X
j=i+1
a
ij
x
(k)
j
, i = 1, 2, ··· , n
(3.43)
为什么说 Gauss Seidel 比 Jacobi 更快呢?从单个元素的递推式可以做一个感性的认知。对于 Jacobi 来说,每一次迭代都是用上
一次的结果,而
Gauss Seidel
则是每一次迭代都用到了最新的结果
9
3.3.3.2 收敛性
同 Jacobi 法推出一个不用求逆的特征多项式
det(λI − G) = det(λI −(D − L)
−1
U)
= det(D − L)
−1
det(λ(D −L) −U)
=
1
det(D −L)
det(λ(D − L) −U)
=
1
det(D −L)
det C
g
(3.44)
于是
推论 3.3.7 (Gauss Seidel 迭代法的特征多项式). 特征多项式等价于 C
g
的行列式为零
C
g
=
λa
11
a
12
··· a
1n
λa
21
λa
22
··· a
2n
.
.
.
.
.
.
.
.
.
.
.
.
λa
n1
λa
n2
··· λa
nn
(3.45)
定理 3.3.8 (由对角占优矩阵判断). 对于方程组 Ax = b,如果
• 若 A 是严格对角占优矩阵,则 Jacobi、Gauss Seidel 迭代法收敛
• 若 A 是弱对角占优矩阵,且为不可约矩阵,则 Jacobi、Gauss Seidel 迭代法收敛
定理 3.3.9 (由正定矩阵判断). 对于方程组 Ax = b
• Jacobi 迭代法收敛 ⇔ A 以及 2D − A 是正定矩阵
• Gauss Seidel 迭代法收敛
(充分条件)
⇐ A 是正定矩阵
这个定理的证明显然不会简单的,但是其中提到了一个一眼瞪出特征多项式的方法值得记住
9
然而,这个所谓的最新的结果我觉得是凑出来的,因为其实也完全可以表示成只和前一次迭代值有关
44
数值分析笔记 第三章 线性方程组解法 Tsui Dik Sang
3.4 超松弛迭代法 (SOR)
3.4.1 方法
把 M 取成
M =
1
ω
(D −ωL) (3.46)
即可,从略
3.4.2 收敛性
引理 3.4.1 (SOR 收敛的必要条件).
0 < ω < 2 (3.47)
定理 3.4.2 (SOR 收敛的充要条件). 对于 Ax = b ,如果
• A 为正定矩阵
• 0 < ω < 2
则 SOR 迭代法收敛
推论 3.4.3 (SOR 收敛的另一个充要条件). 对于 Ax = b ,如果
• A 为严格对角占优矩阵
• 0 < ω ⩽ 1
再来一个超纲的结论
定理 3.4.4 (SOR 收敛速度最快因子). 这个因此应该应该使得
ρ(L
ω
opt
) = min
0<ω<2
ρ(L
ω
) (3.48)
其中 L
ω
= (D − ωL)
−1
[(1 − ω)D + ωU], 具体参加书本 P194 的定义
10
10
再后面的关于块迭代法的内容是真受不了,先不写了。
45
数值分析笔记 第三章 线性方程组解法 Tsui Dik Sang
46
第四章 非线性方程数值解法
4.1 基础定义和定理
这一部分其实都是之前学过的内容了,再提一下。
首先给出零点的定义 (复变里面其实也已经提到过了这个定义了)
定义 4.1.1 (根/零点). 如果有 x
∗
满足线性方程组 f (x
∗
) = 0 那么称 x
∗
为其根
定义 4.1.2 (m 重根). 如果其中 f (x) 可以写成
f(x) = (x − x
∗
)
m
g(x), g(x
∗
) ̸= 0 (4.1)
那么称 x
∗
为方程的 m 重根
判断在一个区间中是否有根,用到的也是高数中出现过的定理
定理 4.1.1 (介值定理). 如果函数 f(x) 在区间 [a, b] 上是连续的,并且 f(a) · f(b) < 0(即 f(a) 和 f (b) 异号),那么在该区
间内至少存在一个零点(即 f (x) = 0 的解)。
4.2 二分法
其基本思想之前也已经接触过了,实际上有一个更通俗的名字:夹逼法。反复二分去夹,最后得到一个较精确的近似解。
4.2.1 近似程度
如果是已知轮数求误差精度,则有
定理 4.2.1 (二分法的误差估计). 设 f(x) 在隔间 [a, b] 上是连续的,且 f(a)f (b) < 0, 则由二分法产生的序列 {x
k
}
+∞
k=0
收敛
于方程 f (x) = 0 的根 x
∗
, 并且误差估计为
|x
k
− x
∗
| ⩽
1
2
k
(b − a) (4.2)
反过来,如果是规定精度,要求总共要进行几轮算法,则有推论
47
数值分析笔记 第四章 非线性方程数值解法 Tsui Dik Sang
推论 4.2.2 (规定精度的对分次数). 若预先给定了根的绝对误差限 ε,要求 |x
k
− x
∗
| ⩽ ε, 则可求得至少要对分的次数是
k ⩾
ln(b − a) − ln ε
ln 2
(4.3)
下面的算法6是给定精度值的算法:同理有给定轮数的算法。
算法 6 二分法
输入: 端点 a,b;
根的误差限 ε;
输出: 近似解 c
1: 由式4.3计算出迭代次数
2: for all n = 1, 2, ··· , k do
3: c =
a+b
2
, 并计算 f (c)
4: if f(c)f (b) < 0 then
5: a = c;
6: else
7: b;
8: end if
9: end for
4.2.2 优点和缺点
4.2.2.1 优点
对 f (x) 的光滑度要求不高,只要其连续即可保证算法收敛。
4.2.2.2 缺点
1. 只能求奇数重实根
1
2. 收敛速度与
1
2
为公比的等比数列相同,还是算慢的!
一般会用其他方法为二分法提供初值。
4.3 不动点法
4.3.1 算法简述
这个方法 (算法7) 就是高中数学数列压轴题里面用到的一个方法: 不动点法:
x
k+1
= φ(x
k
), ( k = 0, 1, 2, ···) (4.4)
最终的解应满足
x
∗
= φ(x
∗
) (4.5)
1
因为只有这样的根在区间两端才能分立于 x 轴上下,对于偶数重根,举一个例子:f(x) = x
2
− 3,如果区间在 [-10,10],就无法使用二分法了
48
数值分析笔记 第四章 非线性方程数值解法 Tsui Dik Sang
算法 7 不动点迭代法
输入: 初始值 x
0
, 根的误差限 ε
输出: 近似解 x
1
1: 将原方程写成 x = φ(x) 的形式
2: while |x
1
− x
0
| < ε do
3: x
1
= φ(x
0
)
4: x
0
← x
1
5: end while
4.3.2 性质
在高中的时候我们是用画图这种的不太严谨的方法来展现该方法收敛求解的过程,下面将给出其收敛到唯一解的充分条件
定理 4.3.1 (简单迭代法唯一收敛的充分条件). 若函数 φ(x) 在 [a, b] 上具有一阶连续导函数,并且满足条件
1. ∀x ∈ [a, b], φ(x) ∈ [a, b]
2.
∃L ∈ (0, 1)
∀x ∈ [a, b]
, |φ
′
(x)| ⩽ Lh,或者写成差分形式 |φ(x
1
) − φ(x
2
)| ⩽ L|x
1
− x
2
|
那么有下面的结论:
推论 4.3.2 (唯一性). 方程 x = φ(x) 在 [a, b] 上有唯一的实根 x
∗
推论 4.3.3 (收敛性). ∀x
0
∈ [a, b],迭代公式4.4收敛,即 lim
k→∞
x
k
= x
∗
推论
4.3.4 (
迭代公式的误差估计
).
|x
k
− x
∗
| ⩽
L
1−L
|x
k
− x
k−1
|
|x
k
− x
∗
| ⩽
L
k
1−L
|x
1
− x
0
|
(4.6)
推论 4.3.5. [迭代函数在解处的导数]
lim
k→∞
x
k+1
− x
∗
x
k
− x
∗
= φ
′
(x
∗
) (4.7)
下面简单说一下证明思路
4.3.2.1 唯一性
构造 g(x) = x −φ( x) ,
• 由条件 1+ 介值定理证得其在 [a, b] 上至少有一实根;
• 由条件 2→ 单调性,证得存在唯一实根 x
∗
49
数值分析笔记 第四章 非线性方程数值解法 Tsui Dik Sang
4.3.2.2 收敛性
4.3.1已经是一个收敛条件了,下面据此再给几个常用的。4.4、4.5 两式相减再结合微分中值定理有
|x
k+1
− x
∗
| = |φ(x
k
) − φ(x
∗
)| = φ
′
(ξ
k
)|x
k
− x
∗
| ⩽ L|x
k
− x
∗
| (4.8)
其中 ξ
k
在 x
k
与 x
∗
之间,然后反复地推,最终得到
0 ⩽ L|x
k+1
− x
∗
| ⩽ ··· ⩽ L
k+1
|x
0
− x
∗
| (4.9)
当 k → ∞ 时原不等式得证。
有时候限定在一个区间内的收敛性比较苛刻,于是就有了在精确解一个小邻域的收敛条件,
定义 4.3.1 (局部收敛性). 设 φ(x) 有不动点 x
∗
,如果
∃R : |x −x
∗
| ⩽ δ, ∀x
0
∈ R, lim
k→∞
x
k
= x
∗
, (4.10)
那么称该迭代法局部收敛的
由于所谓的邻域是人为规定的,所以如果我们选取的区域就是局部收敛要求的邻域,那么其实局部收敛对于我们要研究的问题来
说就是全局收敛的,所以了解一下局部收敛的条件吧!
推论 4.3.6 (局部收敛性). 设 φ(x) 在 x
∗
的一个小邻域内具有一阶连续导函数,并且 |φ
′
(x)| < 1,则不动点法局部收敛。
4.3.2.3 误差估计公式
将公式 4.8 拆成递推公式的形式
•
|x
k
− x
∗
| = |x
k
− x
k+1
+ x
k+1
− x
∗
| ⩽ |x
k
− x
k+1
| + |x
k+1
− x
∗
| ⩽ |x
k+1
− x
k
| + L|x
k
− x
∗
| (4.11)
⇒ |x
k
− x
∗
| ⩽
1
1 − L
|x
k+1
− x
k
| (4.12)
• 对于第二个不等式,只需由微分中值定理证明
|x
k+1
− x
k
| = |φ(x
k
) − φ(x
k−1
)| = |φ
′
(η
k
)(x
k
− x
k−1
)| ⩽ L|x
k
− x
k−1
| (4.13)
然后反复递推,得
|x
k+1
− x
k
| ⩽ ··· ⩽ L
k−1
|x
1
− x
0
| (4.14)
然后结合4.12得证。
4.3.2.4 迭代函数在解处的导数
由式4.8得
x
k+1
− x
∗
x
k
− x
∗
= φ
′
(ξ
k
) (4.15)
并且由 φ
′
(x) 连续,取极限得证原式。
4.3.3 迭代法的收敛阶
二分法的收敛速度显而易见,然而简单迭代法却不是那么容易看出来,不过我们可以类似的进行定义
4.3.3.1 收敛阶的定义和一些判定
50
数值分析笔记 第四章 非线性方程数值解法 Tsui Dik Sang
定义 4.3.2 (收敛阶). 设迭代序列 {x
k
}
+∞
k=0
收敛到 x
∗
,记 e
k
= x
k
− x
∗
, 如果存在常数 c > 0 以及实数 p ⩾ 1,使得
c = lim
k→∞
|e
k+1
|
|e
k
|
p
(4.16)
那么称序列 {x
k
}
+∞
k=0
是 p 阶收敛的
定义 4.3.3 (线性收敛). 当 p=1 时,序列 {x
k
}
+∞
k=0
为线性收敛
推论 4.3.7 (线性收敛的性质). 此时必然有 0 < c < 1
定义 4.3.4 (超线性收敛). 当 p > 1 时,序列 {x
k
}
+∞
k=0
为超线性收敛
推论 4.3.8 (超线性收敛的性质). 此时 p 越大,序列 {x
k
}
+∞
k=0
收敛到 x
∗
的速度越快,
推论 4.3.9 (渐进常数). 两个算法的收敛阶相同时,渐进常数 c 越小,序列 {x
k
}
+∞
k=0
收敛越快
推论 4.3.10 (迭代函数在解处的导数).
lim
k→∞
|x
k+1
− x
∗
|
|x
k
− x
∗
|
= lim
k→∞
|e
k+1
|
|e
k
|
= lim
k→∞
|φ
′
(x
∗
)||e
k
|
|e
k
|
p
= |φ
′
(x
∗
)| (4.17)
4.3.3.2 两种算法的收敛阶讨论
• 显然,二分法的收敛阶就是 1,并且渐进常数 c =
1
2
• 由推论4.3.10可以知道
2
,当 φ
′
(x
∗
) ̸= 0 时,简单迭代法也是线性收敛,并且渐近常数 c = |φ
′
(x
∗
)|
因此,当 |φ
′
(x
∗
)| <
1
2
时简单收敛法就比线性收敛法更优了。
4.3.4 优化的迭代法 *
4.3.4.1 线性收敛序列的 Aitken 加速法
对于充分大的 k,
x
k+1
− x
∗
x
k
− x
∗
≈
x
k+2
− x
∗
x
k+1
− x
∗
≈ c (4.18)
解得
x
∗
≈ x
k
−
(x
k+1
− x
k
)
2
x
k+2
− 2x
k+1
+ x
k
(4.19)
当 k → ∞ 时候上式就取等了,因此就得到了加速公式
2
这就是为什么有这个关于迭代函数在解处的导数的推论的原因
51
数值分析笔记 第四章 非线性方程数值解法 Tsui Dik Sang
定理 4.3.11 (Aitken 加速公式). 定义序列 {y
k
}
+∞
k=0
, 满足
y
k
= x
k
−
(x
k+1
− x
k
)
2
x
k+2
− 2x
k+1
+ x
k
(4.20)
其仍然收敛于 x
∗
, 且收敛速度更快
具体为什么更快的原因超纲了。
4.3.4.2 Steffensen 迭代法
这相当于是对 Aitken 加速法的一个具体实现,给出了具体的迭代公式, 将原来的迭代函数 φ(x) 换成如下函数
ψ(x) = x =
[φ(x) − x]
2
φ(φ(x)) − 2φ(x) + x
(4.21)
这个算法可以达到二阶收敛。具体原因也是超纲了。
4.4 Newton 迭代法
本质上是对二分法的一种改进,但是在形式上是一种特殊的不动点法,通过做切线加速收敛速度,迭代公式如下
x
k+1
= x
k
−
f(x
k
)
f
′
(x
k
)
(4.22)
对于其收敛性也可以用不动点的方法来证明
3
定理 4.4.1 (Newton 迭代法的收敛性). Newton 法在 x
∗
的邻近是平方收敛的
关于牛顿法有很多变形, 这里尝试一语点出其精髓,不再给例题
4.4.1 简化 Newton 法:平行弦方法
也就是固定斜率,偷懒,用第一个点的斜率作为之后的斜率
C =
1
f
′
(x
0
)
x
k+1
= x
k
−
f(x
k
)
C
(4.28)
3
对于牛顿法将其转换成不动点法则
φ(x) = x −
f(x)
f
′
(x)
(4.23)
⇒ φ
′
(x) =
f(x)f
′′
(x)
[f
′
(x)]
2
(4.24)
⇒ φ
′′
(x) =
f
′′
(x
∗
)
f
′
(x
∗
)
(4.25)
如果 x
∗
是 f (x) = 0 的单根,那么有
f(x
∗
) = 0, f
′
(x
∗
) ̸= 0, f
′′
(x
∗
) ̸= 0 · · · (4.26)
所以就知 φ
′
(x
∗
) = 0, φ
′′
(x
∗
) ̸= 0, 也就证明了收敛阶为 2,并且由结论4.3.10可以知道
lim
k→∞
|x
k+1
− x
∗
|
|x
k
− x
∗
|
2
=
f
′′
(x
∗
)
f
′
(x
∗
)
(4.27)
52
数值分析笔记 第四章 非线性方程数值解法 Tsui Dik Sang
4.4.2 Newton 下山法
看算法有点类似于 PID 中的不完全微分算法 即在迭代中再与上次值做一个加权叠加
¯x = x
k
−
f(x
k
)
f
′
(x
k
)
x
k+1
= λ
k
¯x + (1 − λ
k
)x
k
, λ
k
∈ (0, 1)
(4.29)
其中 λ 称为下山因子,
4.4.3 重根情况
对于重根会出现 f
′
(x
∗
) = 0 的情况,然而这只是在精确解时候为零,如果在其附近一阶导不为零的话仍然可以用 Newton
法来求解,
4.4.3.1 重根情形分析
设一种形式来描述重根
定义 4.4.1 (重根情形). 若 f(x) 能表示成
f(x) = (x − x
∗
)
m
g(x), g(x
∗
) ̸= 0 , m ⩾ 2 (4.30)
则 x
∗
为方程 f (x) = 0 的 m 重根
推论 4.4.2 (m 重根的性质).
f(x
∗
) = f
′
(x
∗
) = f
′′
(x
∗
) = ··· = f
(m−1)
(x
∗
) = 0 , f
(m)
(x
∗
) ̸= 0 (4.31)
单只要 f
′
(x
k
) ̸= 0 ,由
|φ
′
(x
∗
)| < 1 (4.32)
则仍然可以用 Newton 法来求解,不过
4
φ
′
(x
∗
) = 1 −
1
m
̸= 0 (4.34)
意味着只有线性收敛,如果想要平方收敛的话还需要优化
4.4.3.2 平方收敛的重根算法
不加证明的给出两种,有时间会给出证明
•
x
k+1
= x
k
− m
f(x
k
)
f
′
(x
k
)
(4.35)
•
x
k+1
= x
k
−
f(x
k
)f
′
(x
k
)
[f
′
(x
k
)]
2
− f(x
k
)f
′′
(x
k
)
(4.36)
4
φ(x) = x −
f(x)
f
′
(x)
= x −
(x − x
∗
)
m
g(x)
(x − x
∗
)
m−1
[mg(x) + g
′
(x)(x − x
∗
)]
(4.33)
但是我还没有证明出来,书上也没有展开证明,有时间我会完善的。
53
数值分析笔记 第四章 非线性方程数值解法 Tsui Dik Sang
4.4.4 类差分插值法
整天求导数对于某些奇怪的函数来说还是显得太繁琐了,如果只有差分形式的加减乘除就会好算一点。所以这类类差分的
Newton 法就是在做这件事
4.4.4.1 弦截法
x
k+1
= x
k
−
f(x
k
)
f(x
k
) − f(x
k−1
)
(x
k
− x
k−1
) (4.37)
4.4.4.2 抛物线法
那就和前面两个点有关
ω = f [x
k
, x
k−1
] + f[x
k
, x
k−1
, x
k−2
](x
k
− x
k−1
)
x
k+1
= x
k
−
2f(x
k
)
ω ±
√
ω
2
−4f[x
k
,x
k−1
]f(x
k
−x
k−1
)
(4.38)
具体证明自己去看书本 p229, 方法用到了牛顿插值法的待定系数过程。
4.5 非线性方程组的数值解法
在这一部分,我们需要用矩阵的思维去看待前面提到的一众非线性解法,并且我觉得最好就是都能亲手编程尝试一下批量化
求解,这才是学这门课的精髓。
4.5.1 非线性方程组的矩阵表示
因此,首先要解决的第一个问题就是用矩阵来表示非线性方程组以及其涉及到的一些计算
定义 4.5.1 (非线性方程组的矩阵表示). 设有非线性方程组
f
1
(x
1
, x
2
, ··· , x
n
) = 0
f
2
(x
1
, x
2
, ··· , x
n
) = 0
.
.
.
f
n
(x
1
, x
2
, ··· , x
n
) = 0
(4.39)
可以写成
F (x) = 0 (4.40)
其中
F (x) =
f
1
(x
1
, x
2
, ··· , x
n
)
f
2
(x
1
, x
2
, ··· , x
n
)
.
.
.
f
n
(x
1
, x
2
, ··· , x
n
)
(4.41)
定义 4.5.2 (雅可比矩阵 (矩阵函数的导数)). 设有非线性方程组
F (x) = 0 (4.42)
54
数值分析笔记 第四章 非线性方程数值解法 Tsui Dik Sang
其雅可比矩阵为
F
′
(x) = J(x) =
∂f
1
∂x
1
∂f
1
∂x
2
···
∂f
1
∂x
n
∂f
2
∂x
1
∂f
2
∂x
2
···
∂f
2
∂x
n
.
.
.
.
.
.
.
.
.
.
.
.
∂f
n
∂x
1
∂f
n
∂x
2
···
∂f
n
∂x
n
(4.43)
4.5.2 不动点迭代法
4.5.2.1 基本使用
扩展一下就行了,无需多言。
定义 4.5.3 (不动点迭代法). 设有非线性方程组
F (x) = 0 (4.44)
可以将其写成
x = φ( x) (4.45)
其中 φ(x) 是一个向量函数,那么可以得到迭代格式
x
(k+1)
= φ(x
(k)
) (4.46)
4.5.2.2 误差分析
使用范数来度量误差
定理 4.5.1 (误差估计).
∥x
∗
− x
(k+1)
∥ ⩽
L
1 − L
∥x
(k)
− x
(k−1)
∥ (4.47)
4.5.2.3 收敛性
定理 4.5.2 (不动点解非线性方程组的收敛性). 设 Φ 在定义域内有不动点 (向量)x
∗
, 则若 Φ 的分量函数有连续偏导数并且
ρ(Φ
′
(x
∗
)) < 1 (4.48)
那么不动点序列迭代法收敛于 x
∗
4.5.3 Newton 迭代法
4.5.3.1 基本使用
定理 4.5.3 (Newton 迭代法). 设有非线性方程组
F (x) = 0 (4.49)
那么可以得到迭代格式
x
(k+1)
= x
(k)
− F
′
(x
(k)
)
−1
F (x
(k)
) (4.50)
在编程上一般使用累加的形式实现 F
′
(x
(k)
)
−1
F (x
(k)
) 这一项。
55
数值分析笔记 第四章 非线性方程数值解法 Tsui Dik Sang
4.5.3.2 收敛性
引理 4.5.4 (方程组的 p 阶收敛性). 一样的,变成矩阵形式而已
lim
k→∞
∥x
(k+1)
− x
∗
∥
∥x
(k)
− x
∗
∥
p
= α, ∃α > 0 (4.51)
那么称 x
(k)
是 p 阶收敛的
定理 4.5.5. 若……,且
∥F
′
(x) − F
′
(x
∗
)∥ < L∥x − x
∗
∥, ∃L < 1 (4.52)
那么 Newton 迭代法至少平方收敛于 x
∗
56
第五章 数值积分与数值微分
重点是数值积分,数值微分方法提到的不多
5.1 数值积分的基础概念
5.1.1 思想与引入
用数值法去解积分而不想算出解析解,通俗来讲就是想要暴力强拆,那么以直代曲、有限分割再求和就是基本的思想。
5.1.1.1 两点积分近似
一个最原式的是通过两点去近似,这当然不精确,不过对于我们进一步探索提供了一个思路。
定义
5.1.1 (
梯形公式
).
设
f
(
x
)
在区间
[
a, b
]
上连续,且边界值已知,则积分值可近似为
ˆ
b
a
f(x)dx ≈
b − a
2
[f(a) + f (b)] (5.1)
另一种近似方法是矩形
定义 5.1.2 (矩形公式). 设 f(x) 在区间 [a, b] 上连续,且边界值已知,则积分值可近似为
ˆ
b
a
f(x)dx ≈ (b − a)f
a + b
2
(5.2)
5.1.1.2 多点扩展
如果在 [a, b] 内选取多个节点 x
k
,利用矩形公式的扩展,我们就可以构造出求积公式
定理 5.1.1 (机械求积). 设 f(x) 在区间 [a, b] 上连续, 则积分值可近似为
ˆ
b
a
f(x)dx ≈
n
X
k=0
A
k
f(x
k
) (5.3)
其中 A
k
为常数,x
k
称为求积节点,A
k
称为求积系数, 亦称为伴随节点 x
k
的权。
5.1.2 代数精度
定义 5.1.3. 如果某个求积公式对于次数不超过 m 的多项式 f (x) 都能准确成立,那么称该求积公式具有 m 次代数精度。
显然根据求积公式的线性性可以得到
57
数值分析笔记 第五章 数值积分与数值微分 Tsui Dik Sang
推论 5.1.2 (求积公式具有 m 次代数精度的充要条件). 需要求积公式对 f(x) = x
k
, k = 0, 1, ··· , m 都成立。即
n
X
k=0
A
k
= b − a
n
X
k
=0
A
k
x
k
=
(b
2
− a
2
)
2
.
.
.
n
X
k=0
A
k
x
m
k
=
(b
m
− a
m
)
m − 1
(5.4)
一开始我以为应该从自由度的角度去理解,给了 n 个点就最多只能有 n-1 的精度,但是后面发现并不一定,这个定理要的只是要
满足尽可能多的方程,并不一定意味着当未知数个数大于方程时一定会出现冲突,有这样的风险,所以有有一个下限5.2.3,但是
可以更高,
推论 5.1.3 (梯形法和矩形法的代数精度). 梯形公式和矩形公式都是有一次代数精度的
5.1.3 求积公式的余项
显然由代数精度的定义直接可以有下面的结论
推论 5.1.4 (求积公式的余项).
R[f ] =
ˆ
b
a
f(x)dx −
n
X
k=0
A
k
f(x
k
) = Kf
(m+1)
(η), η ∈ [a, b] (5.5)
K
是一个不依赖于
f
(
x
)
的常数
当 f (x) 的次数小于等于 m 的多项式时,f
(m+1)
(x) = 0,余项为零,求积公式精确成立,这与代数精度的定义一致。
虽然余项可以推,在实际使用中记住余项意义也不大,但是功利来看,考场上不可能现推吧?所以——下面涉及一些重要余
项系数的时候我会标红的,不一定要记住,但一定要知道在书中哪一页有提到过,从而在开卷考试的时候能快速反应到。
5.1.4 收敛性与稳定性
5.1.4.1 收敛性
定义 5.1.4 (收敛性). 设 f(x) 在区间 [a, b] 上连续,且其导数在该区间上有界如果
ˆ
b
a
f(x)dx = lim
n→∞
n
X
k=0
A
k
f(x
k
) (5.6)
即 lim
n→∞
R[f ] = 0,,则称求积公式5.1.1是收敛的,
5.1.4.2 稳定性
定义 5.1.5 (稳定性). 若
˜
f
k
是由于扰动 δ
k
而产生的近似值,
∀ε > 0, ∃δ > 0, |f (x
k
) −
˜
f
k
| ⩽ δ ⇒ |I
n
(f| − I
n
(
˜
f)| ⩽ ε (5.7)
58
数值分析笔记 第五章 数值积分与数值微分 Tsui Dik Sang
则称求积公式5.1.1是稳定的
这个定义比较不直观,不过记住下面的结论就够了
定理 5.1.5 (稳定性的判定). 若求积公式5.1.1中的系数 A
k
是正数,则该求积公式是稳定的。
5.2 插值型求积
所以这是我们立刻就要引入的第一个数值积分方法,因为机械求积的公式实际上和插值法公式极度相识。
5.2.0.1 求积公式
因此通过多项式积分容易得到
定理 5.2.1 (插值型求积公式). 对于机械积分5.1.1中的系数
A
k
=
ˆ
b
a
l
k
(x)dx, k = 0, 1, ··· , n (5.8)
其中 l
k
是 f (x) 插值函数 L
n
(x) 的系数
5.2.0.2 余项公式
余项也是容易求得的
定理 5.2.2 (插值型求积余项公式). 其中余项 R
n
(f) 为
R
n
(f) =
f
(n+1)
(ξ)
(n + 1)!
ω
n+1
(x), ξ ∈ [a, b] (5.9)
其中 ω
n+1
(x) = (x − x
0
)(x − x
1
) ···(x −x
n
)
插值型求积的积分精度没有上限,但是有下限
推论 5.2.3 (插值型求积公式的精度下限). 形如5.1.1的求积公式至少有 n 次代数精度,如果他是插值型的
5.3 Newton-Cotes 公式
5.3.1 引入
无他,随意选取节点的插值型求积是不好程序实现的,但是对积分区域等分取节点就可以很好的进行规范化编程,这也就是
Newton-Cotes 公式:
定义 5.3.1 (Newton-Cotes 公式). 设将积分区间 [a, b] 划分为 n 等分,步长为 h =
b−a
n
,, 选取等距节点 x
k
= a + kh, k =
0, 1, ··· , n,那么 Newton-Cotes 公式为
I
n
= (b − a)
n
X
k=0
C
(n)
k
f(x
k
) (5.10)
59
数值分析笔记 第五章 数值积分与数值微分 Tsui Dik Sang
其中 C
(n)
k
为 Cotes 系数
其系数可以通过插值法求出
1
推论 5.3.1 (Cotes 系数).
C
(n)
k
=
h
b − a
ˆ
n
0
Y
j=0
j̸=k
t − j
k −j
dt =
(−1)
n−k
nk!(n − k)!
ˆ
n
0
Y
j=0
j̸=k
(t − j)dt (5.11)
可以列出部分 Cotes 系数表如下
n C
(n)
0
C
(n)
1
C
(n)
2
C
(n)
3
C
(n)
4
C
(n)
5
C
(n)
6
C
(n)
7
C
(n)
8
1
1
2
1
2
2
1
6
4
6
1
6
3
1
8
3
8
3
8
1
8
4
7
90
32
90
12
90
32
90
7
90
5
19
288
75
288
50
288
50
288
75
288
19
288
6
41
840
216
840
27
840
272
840
27
840
216
840
41
840
7
751
17280
3577
17280
1323
17280
2989
17280
2989
17280
1323
17280
3577
17280
751
17280
8
989
28350
5888
28350
−928
28350
10496
28350
−4540
28350
10496
28350
−928
28350
5888
28350
989
28350
表 5.1: Cotes 系数表
可以看到在 n=8 时,Cotes 系数已经开始出现负数了,也就意味着该求积公式在高阶的情况有不稳定情况,这也是 Newton-
Cotes 公式的一个固有缺点。
5.3.2 余项
5.3.2.1 奇数阶
没有结论,只能由结论5.2.3确定一个下限
5.3.2.2 偶数
书本上证明其精度可以更高,证明从略
定理 5.3.2 (偶数阶的代数精度). 当阶 n 为偶数是,Newton-Cotes 公式至少具有 n + 1 次代数精度
5.3.3 Simpson 公式
5.3.3.1 定义
是 n=2 时候的 Newton-Cotes 公式
1
但是说实话我还没有亲手证明一次
60
数值分析笔记 第五章 数值积分与数值微分 Tsui Dik Sang
定义 5.3.2 (Simpson 求积公式). 对于 [a,b] 区间上的积分
S =
b − a
6
f(a) + 4f
a + b
2
+ f(b)
(5.12)
5.3.3.2 余项
由定理5.3.2可知 Simpson 公式的代数精度为 3,所以
定理 5.3.3 (Simpson 的余项).
R[f ] = −
b − a
180
b − a
2
f
(4)
(η) (5.13)
5.4 复合求积公式
类比到插值里面就是分段插值公式,就是为了抵抗单函数插值或者求积时候出现的龙格现象,对与求积来说,对抗的是不稳
定现象
5.4.1 复合梯形公式
在每一个子区间上都使用梯形公式
定义 5.4.1 (复合梯形公式).
I ≈ T
n
=
h
2
"
f(a) + 2
n−1
X
k=1
f(x
k
) + f(b)
#
(5.14)
5.4.1.1 余项
对于其余项并不能直接套公式5.1.4,因为其余项是分段的,需要求和得到
2
推论 5.4.1 (复合梯形公式的余项).
R
n
(f) = I − T
n
=
n−1
X
k=0
−
h
3
12
f
′′
(η
k
)
, η
k
∈ (x
k
, x
k+1
)
= −
b − a
12
h
2
f
′′
(η), η ∈ [a, b]
(5.17)
5.4.2 复合 Simpson 公式
我觉得这有一点作弊,由于 Simpson 公式需要三个点,所以每一个区间还需要增加一个点,相当于是在分成 n 份的积分上
再二等分,算出 x
k+
1
2
= x
k
+
1
2
h 的值
2
其中还用到了
min
0⩽k⩽n−1
f
′′
(η
k
) ⩽
1
h
2
n−1
k=0
f
′′
(η
k
) ⩽ max
0⩽k⩽n−1
f
′′
(η
k
) (5.15)
所以
f
′′
(η) =
1
n
n−1
k=0
f
′′
(η
k
), η
k
∈ (x
k
, x
k+1
), η ∈ [a, b] (5.16)
61
数值分析笔记 第五章 数值积分与数值微分 Tsui Dik Sang
定义 5.4.2 (复合 Simpson 公式).
I ≈ S
n
=
h
6
"
f(a) + 4
n−1
X
k=0
f(x
k
) + 2
n−1
X
k=1
f(x
k+
1
2
) + f(b)
#
(5.18)
5.4.2.1 余项
同理的方法可以求和得到余项表达式
推论 5.4.2 (复合 Simpson 公式的余项).
R
n
(f) = I − S
n
= −
(b − a)
180
h
2
4
f
(4)
(η), η ∈ [a, b] (5.19)
所以余项系数应是
1
180
·
1
2
4
=
1
180·2
4
5.5 Romberg 求积公式
这个算法是我在复习的时候才终于理解了的算法,还是比较抽象的,关键是要单独理解好主线的递推化算法和支线的外推
加速算法,然后用封装递推的思想就能理清思路了
5.5.1 主线:梯形公式的递推化
与 Simpson 法相似,也是要增加点,不过使用的拟合公式仍是梯形,于是我们从最初的
n = 0
h = b − a
开始,逐步对区间二
分,得到了梯形的递推化算法 这个算法可以被证明
3
误差提高到了 O(h
2
)
算法 8 梯形的递推化算法
输入: 区间 [a,b], 迭代轮数 n
输出: 将区间分成了 k 份的求积结果 T
k
1: T
0
=
h
2
[f(a) + f (b)]
2: for all k = 1 to n do
3: T
2n
=
1
2
T
n
+
h
2
n−1
X
k=0
f(x
k+
1
2
)
4: end for
推论 5.5.1 (梯形递推的余项).
R
n
(f) = I − T
n
= −
(b − a)
12
h
2
f
′′
(η), η ∈ [a, b] (5.20)
即具有 2 阶精度
3
目前还没有理清楚思路
62
数值分析笔记 第五章 数值积分与数值微分 Tsui Dik Sang
5.5.2 支线:外推加速算法
如果对梯形递推法进一步分割 h,误差肯定会接着下降,但是由于误差阶与 h 有关,可以证明误差阶不会下降,那能否在已
有的分割精度 h 的基础上进一步提高误差阶呢?这种情况之前遇到过,之后常微分方程数值解法的时候也将还会遇到,就是想方
设法消去余项中低阶的 h 项,先看一个引理
引理 5.5.2 (余项系数与 h 的无关性). 设 f(x) ∈ C
∞
[a, b],则有
T (h) = I + α
1
h
2
+ α
2
h
4
+ ··· + α
n
h
2n
+ O(h
2n+2
) (5.21)
其中系数 α
k
与 h 无关
推论 5.5.3 (构造四阶精度). 由 T (h) 与 T
h
2
可以构造出拥有四阶精度的求积公式 S(h), 满足
S(h) =
4T
h
2
− T (h)
3
= I + β
1
h
2
+ β
2
h
4
+ O(h
6
) (5.22)
同理可以由 S(h) 推出拥有六阶精度的求积公式
C(h) =
16S
h
2
− S(h)
15
= I + γ
1
h
2
+ γ
2
h
4
+ O(h
6
) (5.23)
于是我们可以归纳加速公式了
定理 5.5.4 (Richardson 外推加速公式).
T
m
(h) =
4
m
4
m
− 1
T
m
−
1
h
2
−
1
4
m
− 1
T
m
−
1
(h) (5.24)
推论 5.5.5 (Richardson 外推加速公式的余项).
T
m
(h) = I + δ
1
h
2(m+1)
+ δ
2
h
2(m+2)
+ O(h
2(m+3)
) (5.25)
而 Romberg 算法就是递推化公式以及、Richardson 外推加速公式的结合
算法 9 Romberg 算法
输入: 区间 [a,b], 迭代轮数 n
输出: 将区间分成了 k 份的求积结果 T
k
1: T
0
=
h
2
[f(a) + f (b)]
2: for all k = 1 to n do
3: 按照算法8计算 T
(k)
0
,
4: end for
5: while |T
(0)
k
− T
(0)
k−1
< ε| do
6: 按照 T
(k)
k+1
=
4
k
T
(k)
k
−T
(k)
k−1
4
k
−1
计算出整个 T 表
7: end while
63
数值分析笔记 第五章 数值积分与数值微分 Tsui Dik Sang
5.6 自适应算法
刚刚的所有二分方法都是无脑二分,对于每一个区间都要细分,这显然是对计算量不友好的,因此如果能对一些可能会偏差
较大 (比如波动变化剧烈) 的区间细分更多,而对平缓区域细分更少,可以大大降低计算量。
这个方法需要“判断”,因此在实操上我自己觉得没有这么容易,上面就是思路,具体的代码在网站上用 ai 写的,但是具体
算法 10 自适应算法
输入
:
区间
[a,b],
精度
ε
输出: 复合误差要求的结果 S
n
, 其中 n 是某一个区间分割的最高阶数
1: n=1,S
0
= I, (I 是可能题目有给的精确值,如果没有的话,下面的 while 从 1 开始)
2: S
1
= S
1
(a, b) =
h
6
[f(a) + 4f
a+b
2
+ f(b)]
3: while
S
n−1
(a, b) −
S
n
(a,
a+b
2
) + S
n
(
a+b
2
, b)
> ε (对于左边区间来说) do
4: 接着分!
5: n++
6: end while
7: 对于右边也是一样!
内容还未看懂
5.7 Gauss 求积公式
重新审视机械求积公式5.1.1,
ˆ
b
a
f(x)dx ≈
n
X
k=0
A
k
f(x
k
) (5.26)
可以发现假如 x
k
的选取也是未知的话,理应有 2n + 2 个未知参数,也就是说理论上求积精度可以达到 2n + 1 次,而前面我们
一直在用等距节点,相当于就抹杀了这些精度,使得最高精度只能到 n。所以 Gauss 方法就通过选取合适的节点 x
k
来把这些精
度发挥出来
5.7.1 二点插值
先从简单的两个点的情况来看,
对于插值公式
ˆ
1
−1
f(x)dx ≈ A
0
f(x
0
) + A
1
f(x
1
) (5.27)
试确定四个未知数 A
0
, A
1
, x
0
, x
1
,使得其具有尽可能的高的代数精度
由结论5.1.2可知,精度越高就要满足方程组中更多的方程,这里有四个未知数,于是可以列四个方程
A
0
+ A
1
= 2
A
0
x
0
+ A
1
x
1
= 0
A
0
x
2
0
+ A
1
x
2
1
=
2
3
A
0
x
3
0
+ A
1
x
3
1
= 0
(5.28)
64
数值分析笔记 第五章 数值积分与数值微分 Tsui Dik Sang
求解得到
A
0
= 1
A
1
= 1
x
0
= −
1
√
3
x
1
=
1
√
3
(5.29)
所以
推论 5.7.1 (二点 Gauss 求积公式具有 3 阶精度).
ˆ
1
−1
f(x)dx ≈ f
−
1
√
3
+ f
1
√
3
(5.30)
的代数精度为 3。
5.7.2 权函数 ρ(x ) 的引入
书中并没有说明白为什么要引入权函数 ρ(x)
4
,
定义 5.7.1 (Gauss 求积公式). 对于求积公式
ˆ
b
a
f(x)ρ(x)dx ≈
n
X
k=0
A
k
f(x
k
) (5.31)
如果其具有 2n + 1 次代数精度,则其节点 x
k
称为 Gauss 节点,相应的公式称为 Gauss 求积公式,
关于如何证明存在性等书中 p118-120 给了一系列推导,不过其实列方程解方程5.1.2已经可求了,下面给一些结论
推论
5.7.2 (
系数的正定性
).
Gauss
求积公式的系数
A
k
都是正数
推论 5.7.3 (Gauss 公式的稳定性). Gauss 求积公式是稳定的
5.7.3 Gauss-Legendre 公式
定义 5.7.2 (Gauss-Legendre 公式). 即当权函数 ρ(x) = 1 时候的求积公式,并且求积范围为 [−1, 1],
ˆ
1
−1
f(x)dx ≈
n
X
k=0
A
k
f(x
k
) (5.32)
下面的结论不做证明。
推论 5.7.4 (Gauss-Legendre 公式的节点). Legendre 多项式 P
n+1
(x) 的零点就是求积公式的高斯点
关于前五次 Gauss-Legendre 公式的节点和系数参见下面的表格 以后做题看到也是直接查表即可。然后一直到 p126 的内容都从
略
4
ai 说是是为了处理一些特殊形式,类似于一种待定系数法了
65
数值分析笔记 第五章 数值积分与数值微分 Tsui Dik Sang
表
5.2:
前五次
Gauss-Legendre
公式的节点和系数
n 节点 x
k
系数 A
k
0 0 2
1 ±0.5773503 1
2
±0.7745967
0
0.5555556
0.8888889
3
±0.8611363
±0.3399810
0.3478548
0.6521452
4
±0.9061798
±0.5384693
0
0.2369269
0.4786287
0.5688889
5
±0.9324699
±0.6612094
±0.2386192
0.1713245
0.3607616
0.4679139
5.8 多重积分
嵌套用前面所学的求积公式即可,还是要多做题.
5.9 数值微分
5.9.1 中点方法
使用数值微分的差分公式直接就可以得到一个数值微分解法
定理 5.9.1.
f
′
(a) ≈
f(a+h)−f (a)
h
f
′
(a) ≈
f(a)−f (a−h)
h
⇒ f
′
(a) ≈
f(a + h) − f( a − h)
2h
(5.33)
前面两个公式的误差阶都是 O(h)。而后一个公式的误差阶是 O(h
2
)
5.9.2 插值型求导公式
与插值型求积是一样的,也就是先拟合再求导。
定义 5.9.1 (插值型求导公式).
f
′
(x) ≈ P
′
n
(x) (5.34)
推论 5.9.2 (插值型求导公式的余项).
R
n
(f) =
f
(n+1)
(ξ)
(n + 1)!
ω
′
n+1
(x) +
ω
n+1
(x)
(n + 1)!
d
dx
f
(n+1)
(ξ), ξ ∈ [a, b] (5.35)
具体方法不再赘述,回归插值法老套路,当然,插值法不止有拉格朗日插值法,还有三次样条……同样可以延伸,但是那这不是
这一章的重点了
66
第六章 常微分方程初值问题数值解法
6.1 引入
首先重温一下高数中对常微分方程初值问题有唯一解的定理
定理 6.1.1 (Lipchitz 条件). 对于一阶常微分方程:
y
′
= f (x, y), x ∈ [a, b]
y(x
0
) = y
0
(6.1)
如果 ∃L > 0, 使得
|f(x, y
1
) − f(x, y
2
)| ⩽ L|y
1
− y
2
|, ∀x ∈ [a, b], y
1
, y
2
∈ R (6.2)
那么在区间 [a, b] 上存在唯一解 y(x), 上面满足的条件就是 Lipschitz 条件
推论 6.1.2 (解关于初值条件的敏感性). 设处置问题
y
′
= f (x, y), x ∈ [a, b]
y(x
0
) = s
(6.3)
的解为 y(x, s), 则
|y(x, s
1
) − y(x, s
2
)| ⩽ e
L|x−x
0
|
|s
1
− s
2
|, ∀x ∈ [a, b], s
1
, s
2
∈ R (6.4)
这里数值解法的思路和求积有一点类似,就是分割点,形式上有两种:单步法和多步法
6.1.1 单步法
6.1.2 Euler 法
6.1.2.1 方法推导
这个方法的思路应该是从6.3来的,先取一个满足方程的点作为迭代初始点
1
下一个要求的点是 x
1
= x
0
+ h, 如果 h 取得比
较小的话,可以用导数求近似差分式,即
y
n+1
− y
n
x
n+1
− x
n
= f (x
n
, y
n
) (6.5)
于是就得到了递推公式
1
一般这是左边界点
67
数值分析笔记 第六章 常微分方程初值问题数值解法 Tsui Dik Sang
定理 6.1.3 (Euler 法).
y
n+1
= y
n
+ h · f(x
n
, y
n
), n = 0, 1, ··· , N − 1 (6.6)
Euler 法也可以从另一个角度推出
直接对6.3两边在 [x
n
, x
n+1
] 上积分,得到
y(x
n
+ 1) = y(x
n
) +
ˆ
x
n+1
x
n
f(t, y(t))dt (6.7)
对积分项做矩形近似的话就是 Euler 法的递推式了 (x
n+1
− x
n
= h)
6.1.2.2 后退的 Euler 法
如果选的是后一个点的导数或者是后一个点做参考做矩形近似就得到了后退的 Euler 法
定理 6.1.4 (后退的 Euler 法).
y
n+1
= y
n
+ hf(x
n+1
, y
n+1
) (6.8)
虽然推导的方法类似,不过后退的 Euler 法是隐式的!书上 [2] 说这在一些情况下更有好处,不过至少我现在还没有看出来
6.1.2.3 误差
Euler 法属于单步法, 也就是说每一个后续都只有前面的一个点的值来决定,因此误差必然是很大的,尤其越到后面的点偏
移越离谱。如果单看第二个点的误差,可以使用泰勒展开,
推论 6.1.5 (Euler 法的误差).
y(x
n+1
) − y
n+1
=
h
2
2
y
′′
(ξ
n
) ≈
h
2
2
f
′′
(x
n
) (6.9)
再后面的误差就没眼看了……
6.1.3 梯形方法
刚刚的欧拉法是用矩形去近似积分推导的,如果用梯形去,应该会更精确,这就是这个方法的改进。
定义 6.1.1 (梯形法).
y
(0)
n+1
= y
n
+ hf(x
n
, y
n
)
y
(k+1)
n+1
= y
n
+
h
2
[f(x
n
, y
n
) + f(x
n+1
, y
(k)
n+1
)], k = 0, 1, 2, ···
(6.10)
算梯形公式一定要用到后面的点,因此梯形法是隐式的。显然,他的误差肯定会比欧拉法小。
6.1.4 改进的 Euler 法
隐式算法的缺点就是要迭代,这会大大提高计算量,如果结合精度不够但不用隐式计算的 Euler 法预测一个初值,在用初值
迭代一步得到迭代终点值,那么这个迭代的精准值精确度应该会和由一个随便设的初值 y
(0)
n+1
迭代好多次的结果差不多,但是计
算量降低了,这就是改进的 Euler 法
68
数值分析笔记 第六章 常微分方程初值问题数值解法 Tsui Dik Sang
定理 6.1.6 (改进的 Euler 法).
y
(0)
n+1
= y
n
+ hf(x
n
, y
n
)
y
n+1
= y
n
+
h
2
[f(x
n
, y
n
) + f(x
n+1
, y
(0)
n+1
)]
(6.11)
或者写成平均值的形式
y
p
= y
n
+ hf(x
n
, y
n
)
y
c
= y
n
+ hf(x
n+1
, y
p
)
y
n+1
=
y
p
+y
c
2
(6.12)
这里做一次迭代即可了
2
6.1.5 单步法总结及误差分析
根据前面的几个方法,于是我们可以先给一个定义
定义
6.1.2 (
单步法
).
求解下一个点数值解的公式如下的方法称为
单步法
,
y
n+1
= y
n
+ hφ(x
n
, y
n
, y
n+1
, h) (6.13)
然后开始分析误差
定义 6.1.3 (局部截断误差).
T
n+1
= y(x
n+1
) − y(x
n
) − hφ(x
n
, y
n
, y
n+1
, h) (6.14)
“局部”是因为要假设在 x
n
点以及之前是没有误差的,上式的结果才是严格的误差。
因此可以粗略理解为每一步额外产生的误差。
定义 6.1.4 (P 阶精度与误差主项). 将局部截断误差可以写成下面的形式:
T
n+1
= O(h
p+1
) = ψ(x
n
, y(x
n
))h
p+1
+ O(h
p+2
), P ⩾ 1 (6.15)
其中 P 称为精度阶数,O(h
p+1
) 称为误差主项
一般使用泰勒展开可以得到精度的阶数
推论
6.1.7 (
Euler
法的精度
). Euler
法的精度是
1
阶的。
推论 6.1.8 (改进的 Euler 法、梯形法的精度). 改进的 Euler 法、梯形法的精度是 2 阶的。
这是一般的方法,如果碰到一些奇奇怪怪的求解公式要让其阶数尽可能的高,只需要让其低次 h 项尽可能系数为零就行。
比如 14-15 东南大学的 18 题试确定系数 A, B 使得求解公式的误差阶数尽可能的高
y
n+1
= y
n
+ h
Af(x
i
, y
i
) + (1 − A)f(x
i
+ Bh, y
i
+
4
5
hf(x
i
, y
i
))
(6.16)
2
书上将 y
(0)
n+1
写做 ¯y
n+1
69
数值分析笔记 第六章 常微分方程初值问题数值解法 Tsui Dik Sang
首先 Taylor 展开所有能展开的项
R
i+1
= y(x
i+1
) − y(x
i
) − h
Af(x
i
, y
i
) + (1 − A)f(x
i
+ Bh, y
i
+
4
5
hf(x
i
, y
i
))
= y(x
i
) + hy
′
(x
i
) +
h
2
2
y
′′
(x
i
) +
h
3
6
y
′′′
(x
i
) + O(h
4
) − hAy
′
(x
i
)
− h(1 − A)
f(x
i
, y(x
i
)) +
∂f(x
i
, y(x
i
))
∂x
Bh +
∂f(x
i
, y(x
i
))
∂y
4
5
hy
′
(x
i
)
+
1
2
∂
2
f(x
i
, y(x
i
))
∂x
2
B
2
h
2
+
1
2
∂
2
f(x
i
, y(x
i
))
∂y
2
4
5
h
2
y
′
(x
i
)
2
+
∂
2
f(x
i
, y(x
i
))
∂x∂y
4
5
hBhy
′
(x
i
) + O(h
3
)
#
= h · 0 + h
2
1
2
− (1 − A)B
∂f(x
i
, y(x
i
))
∂x
+
1
2
−
4
5
(1 − A)
∂f(x
i
, y(x
i
))
∂
y
′
i
(x
i
)
+ h
3
"
1
6
y
′′′
(x
i
) −
1
2
(1 − A)B
2
∂
2
f(x
i
, y(x
i
))
∂x
2
−
1
2
4
5
2
(1 − A)
∂
2
f(x
i
, y(x
i
))
∂y
2
y
′
(x
i
)
2
−
4
5
(1 − A)
∂
2
f(x
i
, y(x
i
))
∂x∂y
Bhy
′
(x
i
)
#
(6.17)
于是要想阶数尽可能高
(a − A)B =
1
2
4
5
(1 − A) =
1
2
(6.18)
解得
A =
3
8
B =
4
5
(6.19)
此时具有三阶精度
6.2 Runge-Kutta 方法
6.2.1 一般形式
多步法肯定会比单步法误差小。然而这里需要声明一点:Runge-Kutta 法不是多步法!具体原因我会在介绍线性多步法的时
候详细说明
3
。那么我们从哪里插入节点呢?抓手还是 Euler 法,注意到式6.7中,我们直接就将积分项写成关于一个点的函数,
如果想要增加精确度,我们可以在这里做文章,插入更多的点,结合上一章数值求积分的知识我们可以找到插点公式:
ˆ
x
n+1
x
n
f(t, y(t))dt ≈ h
n
X
i=1
c
i
f(x
n
+ λ
i
h, y
n
+ µ
i
h) (6.20)
于是 Runge-Kutta 法就来了
定理 6.2.1 (r 级显式 Runge-Kutta 法).
y
n+1
= y
n
+ hφ(x
n
, y
n
, y
n+1
, h) (6.21)
其中
φ(x
n
, y
n
, y
n+1
, h) =
s
X
i=1
c
i
K
i
K
1
= f (x
n
, y
n
)
K
i
= f
x
n
+ λ
i
h, y
n
+ h
i−1
X
j=1
µ
ij
K
j
, i = 1, 2, ··· , s
(6.22)
显然,r=1 的时候的一阶 Runge-Kutta 法就是 Euler 法,
3
然而在还没有理解透后面的线性多步法之前,我一直以为 Runge-Kutta 法是多步法!
70
数值分析笔记 第六章 常微分方程初值问题数值解法 Tsui Dik Sang
6.2.2 二阶显式 Runge-Kutta 法
6.2.2.1 基本定义
也就是说要插多入一个点,
定义 6.2.1 (二阶显式 Runge-Kutta 法).
y
n+1
= y
n
+ h(c
1
K
1
+ c
2
K
2
)
K
1
= f (x
n
, y
n
)
K
2
= f (x
n
+ λ
2
h, y
n
+ µ
21
hK
1
)
(6.23)
其中 c
1
, c
2
, λ
2
, µ
21
均为待定常数,合适选取这些系数可以使得阶数 p 尽可能的高
6.2.2.2 常数的确定
首先来看截断误差
T
n+1
= y(x
n+1
) + y(x
n
) − h[c
1
f(x
n
, y
n
) + c
2
f(x
n
+ λ
2
h, y
n
+ µ
12
hf
n
)] (6.24)
同样使用泰勒展开来化简式子,注意的是这里需要使用二元泰勒展开
4
代入一波化简有
T
n+1
= (1 − c
1
− c
2
)f
n
h +
1
2
− c
2
λ
2
f
′
x
(x
n
, y
n
)h
2
+
1
2
− c
2
µ
21
f
′
y
(x
n
, y
n
)f
n
h
2
+ O(h
3
) (6.27)
然后我们去消阶即可,把二阶都消掉,需要
1 − c
1
− c
2
= 0
1
2
− c
2
λ
2
= 0
1
2
− c
2
µ
21
= 0
(6.28)
显然,这个解不唯一,不过可以写成关于 c
2
的形式
c
1
= 1 − c
2
λ
2
=
1
2c
2
µ
21
=
1
2c
2
c
2
̸= 0 (6.29)
定义 6.2.2 (中点公式). 如果 c
2
= a, 得到的就是中点公式
y
n+1
= y
n
+ hK
2
K
1
= f (x
n
, y
n
)
K
2
= f
x
n
+
h
2
, y
n
+
h
2
K
1
(6.30)
4
y(x
n+1
= y
n
+ hy
′
n
+
h
2
2
y
′′
n
+
h
3
3!
y
′′
n
+ O(h
4
)) (6.25)
y
′
n
= f (x
n
, y
n
) = f
n
y
′′
n
=
d
dx
f(x
n
, y(x
n
)) = f
′
x
(x
n
, y
n
) + f
′
y
(x
n
, y
n
)f
n
y
′′′
n
=
d
2
dx
2
f(x
n
, y(x
n
)) = f
′′
x
(x
n
, y
n
) + f
′′
y
(x
n
, y
n
)f
n
+ f
′
y
(x
n
, y
n
)f
′
x
(x
n
, y
n
) + f
′
y
(x
n
, y
n
)f
′
y
(x
n
, y
n
)f
n
f(x
n
+ λ
2
h, y
n
+ µ
21
hf
n
) = f (x
n
, y
n
) + λ
2
hf
′
x
(x
n
, y
n
) + µ
21
hf
′
y
(x
n
, y
n
)f
n
+
(λ
2
h)
2
2
f
′′
x
(x
n
, y
n
) +
(µ
21
h)
2
2
f
′′
y
(x
n
, y
n
)f
2
n
+ O(h
3
)
(6.26)
71
数值分析笔记 第六章 常微分方程初值问题数值解法 Tsui Dik Sang
6.2.3 三阶、四阶显式 Runge-Kutta 法
6.2.3.1 三阶显式 Runge-Kutta 法
再增加一个点
定义 6.2.3 (三阶显式 Runge-Kutta 法).
y
n+1
= y
n
+ h(c
1
K
1
+ c
2
K
2
+ c
3
K
3
)
K
1
= f (x
n
, y
n
)
K
2
= f (x
n
+ λ
2
h, y
n
+ µ
21
hK
1
)
K
3
= f (x
n
+ λ
3
h, y
n
+ µ
31
hK
1
+ µ
32
hK
2
)
(6.31)
其中 c
1
, c
2
, c
3
, λ
2
, λ
3
, µ
21
, µ
31
, µ
32
均为待定常数,
其最高精度的参数选取证明类似的:
推论 6.2.2 (三阶显式 Runge-Kutta 法的 4 阶精度情况).
c
1
+ c
2
+ c
3
= 1
λ
2
= µ
21
λ
3
= µ
31
+ µ
32
c
2
λ
2
+ c
3
λ
3
=
1
2
c
2
λ
2
2
+ c
3
λ
2
3
=
1
3
c
3
λ
2
µ
32
=
1
6
(6.32)
同样有很多解,并且表示不完,这里给一个常见的公式
推论 6.2.3 (一个常见的三阶显式 Runge-Kutta 法).
y
n+1
= y
n
+
h
6
(K
1
+ 4K
2
+ K
3
)
K
1
= f (x
n
, y
n
)
K
2
= f
x
n
+
h
2
, y
n
+
h
2
K
1
K
3
= f (x
n
+ h, y
n
− hK
1
+ 2hK
2
)
(6.33)
6.2.3.2 四阶显式 Runge-Kutta 法
72
数值分析笔记 第六章 常微分方程初值问题数值解法 Tsui Dik Sang
定义 6.2.4 (四阶显式 Runge-Kutta 法).
y
n+1
= y
n
+ h(c
1
K
1
+ c
2
K
2
+ c
3
K
3
+ c
4
K
4
)
K
1
= f (x
n
, y
n
)
K
2
= f (x
n
+ λ
2
h, y
n
+ µ
21
hK
1
)
K
3
= f (x
n
+ λ
3
h, y
n
+ µ
31
hK
1
+ µ
32
hK
2
)
K
4
= f (x
n
+ λ
4
h, y
n
+ µ
41
hK
1
+ µ
42
hK
2
+ µ
43
hK
3
)
(6.34)
其中 c
i
, λ
i
, µ
ij
均为待定常数,
参数选取方程更加复杂了,这里给出一个常用的四阶显式 Runge-Kutta 法
推论 6.2.4 (四阶显式 Runge-Kutta 法).
y
n+1
= y
n
+
h
6
(K
1
+ 2K
2
+ 2K
3
+ K
4
)
K
1
= f (x
n
, y
n
)
K
2
= f
x
n
+
h
2
, y
n
+
h
2
K
1
K
3
= f
x
n
+
h
2
, y
n
+
h
2
K
2
K
4
= f (x
n
+ h, y
n
+ hK
3
)
(6.35)
6.2.4 变步长 Runge-Kutta 法
我一开始以为是要自适应,然而实际上书中并没有搞这么深,仅简单分析了一下二分步长对精度的影响。我们以四阶 Runge-
Kutta 法为例,,公式的局部截断误差为 O(h
5
),所以
y(x
n+1
) − y
(
h
2
)
n+1
≈ ch
5
(6.36)
如果折半步长,即有两步误差,每一步都是 c
h
2
5
, 所以
推论 6.2.5 (四阶 Runge-Kutta 法二分后的精度).
y(x
n+1
) − y
(
h
2
)
n
+1
≈ 2c
h
2
5
=
h
5
16
(6.37)
也就是说步长折半后误差减小到了原来的
1
16
,
6.3 单步法误差分析
6.3.1 收敛性
6.3.1.1 基础定义
先给一个抽象的一般性定义
73
数值分析笔记 第六章 常微分方程初值问题数值解法 Tsui Dik Sang
定义 6.3.1 (收敛性). 若一个但不算法对于固定的 x
n
= x
0
+ nh,当 h → 0 时候又 y
n
→ y(x
n
) 即数值解趋近于准确解,那
么这个算法称为收敛的。
显然,局部截断误差以及由此定义的 p 阶精度不足以证明收敛性,
定义 6.3.2 (整体截断误差).
R
n
(y
n
) = y(x
n
) − y
n
(6.38)
其中 T
k+1
是局部截断误差,y(x
n
) 是该点的准确解,y
n
是数值解。
显然,如果整体误差 ∼ O(h
p
),那么该数值方法收敛是肯定的,下面有一个定理:
定理 6.3.1 (整体截断误差精度定理). 假设单步法 (6.1.2)
• 具有 p 阶精度,
• 且增量函数 φ(x, y, h) 关于 y 满足 Lipschitz 条件,即
|φ(x, y, h) − φ(x, ¯y, h)| ⩽ |y − ¯y| (6.39)
• 初值 y
0
是准确的,即 y
0
= y(x
0
)
则其整体截断误差满足
y(x
n
) − y
n
= O(h
p
) (6.40)
证明看起来好烦的样子,所以从略。
推论 6.3.2 (单步法收敛条件).
定理6.3.1中的p ⩾ 1 ⇔ 该单步法是收敛的 (6.41)
6.3.2 相容性
这是一个由收敛性引申出来的概念,看定义就知道。
定义 6.3.3 (相容性). 若单步法 (6.1.2) 的增量函数满足
φ(x, y, 0) = f(x, y) (6.42)
或者换一种形式
lim
h→0
y
1
= f (x, y) (6.43)
显然,第二种形式就容易理解了,容易知道相容性与收敛性等价,因此6.3.2也可以对收敛性成立
推论 6.3.3 (单步法收敛条件).
定理6.3.1中的p ⩾ 1 ⇔ 该单步法是相容的 (6.44)
74
数值分析笔记 第六章 常微分方程初值问题数值解法 Tsui Dik Sang
6.3.3 稳定性
在前面多种数值方法中,不止一次提到过这个词,最通用的定义是定义1.2.1 我们可以由此给出对于常微分方程初值问题求
解算法的算法稳定性定义
定义 6.3.4 (常微分方程初值问题算法的稳定性). 若一个数值方法在节点值 y
0
上存在 δ 的扰动,于以后各节点值 y
m
(m > n)
上产生的偏差均不超过 δ, 则称这个算法是稳定的
6.3.3.1 模型方程 y
′
= λy 的建立
首先要证明之前给出的方程形式 y
′
= f (x, y) 均能转化成或者在“稍微”近似后能够转化成 y
′
= λy 的形式,其中 λ 是复数。
• 这个在邻域 (¯x, ¯y) 展开之后略去高阶项
5
y
′
= f (x, y) = f(¯x, ¯y) + f
′
x
(x, y)(x − ¯x) + f
′
y
(x, y)(y − ¯y) + O((x − ¯x)
2
+ (y − ¯y)
2
) (6.46)
然后再稍微做线性变换就可以得 u
′
= λu 的形式,虽然变了因变量形式,但是不影响通过这个 λ 来限定稳定区间。
• 对于微分方程组 y
′
= Ay
6
, 我们将在刚性方程那里讲,形式上有相似之处
既然有了线性的模型方程,那么对于不同方法就都好办了。
6.3.3.2 Euler 法的稳定性
对于模型方程的递推公式如下
y
n+1
= (1 + hλ)y
n
(6.47)
如此一来,误差传递函数也是线性的了,
ε
n+1
= (1 + hλ)ε
n
(6.48)
如果要求满足稳定性条件 |y
n+1
| ⩽ |y
n
|,则
|1 + hλ| ⩽ 1 (6.49)
⇒ −2 < hλ < 0 (6.50)
这在复平面上是一个圆域,在 λ 固定的情况下,合适选取步长 h 使得 hλ 落在这个区域上,那么 Euler 法就是稳定的。
定义 6.3.5 (绝对稳定域). 单步法解模型方程若得到
y
n+1
= E(hλ)y
n
(6.51)
满足 |E(hλ)| ⩽ 1, 则称 E(hλ) 的值域为绝对稳定域其与实轴的交点称为绝对稳定区间,
5
这个“略去”十分玩味啊,后面我笔者发现就是指线性部分,比如对于
y
′
= −100(y − x
2
) + 2x (6.45)
中的 λ 就直接取 −100 了,
6
这里 A 是 m × m 的 Jacobi 矩阵
75
数值分析笔记 第六章 常微分方程初值问题数值解法 Tsui Dik Sang
6.3.3.3 二阶 Runge-Kutta 法的稳定性
这里选取的是中点公式时候的二阶 Runge-Kutta 法 (方程6.30),同样可以得到
y
n+1
=
1 + hλ +
(hλ)
2
2
y
n
(6.52)
于是
E(λh) = 1 + hλ +
(hλ)
2
2
(6.53)
代入绝对稳定条件 |E(hλ)| ⩽ 1, 解得
−2 < hλ < 0 (6.54)
6.3.3.4 更高阶的 Runge-Kutta 法的稳定性
类似的对于三阶四阶
E
3
(hλ) = 1 + hλ +
(hλ)
2
2!
+
(hλ)
3
3!
E
4
(hλ) = 1 + hλ +
(hλ)
2
2!
+
(hλ)
3
3!
+
(hλ)
4
4!
(6.55)
解得
三阶: − 2.51 < hλ < 0
四阶: − 2.78 < hλ < 0
(6.56)
6.3.3.5 隐式 Euler 法的稳定性
以这个简单的特例来考察隐式法的稳定性,对于隐式
Euler
法,由模型方程得到的递推公式应该是
y
n+1
= y
n
+ λhy
n+1
⇒ y
n+1
=
y
n
1 − hλ
(6.57)
所以
E(hλ) =
1
1 − hλ
(6.58)
可得其稳定域是 (1,0) 为圆心的单位圆外部,稳定域是
h ∈ (0 , +∞), λ < 0 (6.59)
6.3.3.6 梯形法的稳定性
可得此时
E(hλ) =
1 +
hλ
2
1 −
hλ
2
(6.60)
若满足条件发现
−∞ < hλ < 0 (6.61)
如果 λ < 0 稳定区间 h ∈ (0, +∞)
7
这也就是意味着隐式 Euler 法和梯形法是绝对稳定的,在选取步长时无需考虑稳定性问题,
只需考虑精度等,我们定义其性质
定义 6.3.6 (A-稳定). 若一个数值方法的稳定域包含整个左半平面 {hλ|Re(hλ) < 0},则称这个方法是 A-稳定的
7
很神奇,居然稳定域这么大,但是我觉得本质上是因为梯形法在卡 bug,他是隐式的,实际上只能由迭代法求 y
n+1
,如果 y
n+1
是绝对精准的,那么当然
h 怎么样取都是稳定的,然而实际上迭代法决定了 y
n+1
的值是有误差的,因此这个一直稳定的情况也就只出现在理想区间里面
76
数值分析笔记 第六章 常微分方程初值问题数值解法 Tsui Dik Sang
6.4 线性多步法
6.4.1 一般公式与基础概念
直接上定义,就能知道为什么 Runge-Kutta 法不是多步法了
定义 6.4.1 (线性多步法). 如果计算 y
n+k
时,除用 y
n+k−1
的值,还用到了 Y
n+i
(i = 0, 1, ··· , k −2) 的值,那么这个方法称
为线性多步法,其一般公式为
y
n+k
=
k−1
X
i=0
α
i
y
n+i
+ h
k−1
X
i=0
β
i
f
n+i
(6.62)
其中
• y
n+i
为 y(x
n+1
) 的近似
• f
n+i
= f (x
n+i
, y
n+i
)
• x
n+1
= x
n
+ ih
推论 6.4.1 (显式 k 步法和隐式 k 步法). • 如果 β
k
= 0, 则称6.62为显式 k 步法,此时 Y
n+k
可以直接由递推公式得出
• 如果 β
k
̸= 0, 则称6.62为隐式 k 步法,需要用到迭代法来求解 y
n+k
线性多部分也有局部误差,定义思路是一样的,但是形式还是有一点不同
定义 6.4.2 (局部截断误差).
T
n+k
= y(x
n+k
) −
k−1
X
i=0
α
i
y
n+i
− h
k−1
X
i=0
β
i
y
′
(x
n+i
) (6.63)
上面的定义中我们通过微分方程给的条件将 f
n+i
替换成了 y
′
(x
n+i
),
定义归定义,我们还是通过泰勒展开来处理成高阶线性和的形式
y(x
n
+ ih) = y(x
n
) + ihy
′
(x
n
) +
(ih)
2
2
y
′′
(x
n
) + O(h
3
)
y
′
(x
n
+ ih) = y
′
(x
n
) + ihy
′′
(x
n
) + O(h
2
)
(6.64)
代入得到
推论 6.4.2 (局部误差的幂级表示).
T
n+k
= c
0
y(x
n
) + c
1
hy
′
(x
n
) + c
2
h
2
y
′′
(x
n
) + ··· + c
k−1
h
k−1
y
(k)
(x
n
) + O(h
k
) (6.65)
其中由待定系数有
c
0
= 1 − (α
0
+ α
1
+ ··· + α
k−1
)
c
1
= k − [α
1
+ 2α
2
+ ··· + (k − 1)α
k−1
] − (β
0
+ β
1
+ ··· + β
k−1
)
.
.
.c
q
=
1
q !
[k
p
− (α
1
+ 2
q
α
2
+ ··· + (k − 1)
q
α
k−1
)] −
1
(q−1)!
[β
1
+ 2
q − 1
β
2
+ ··· + k
q − 1
β
k−1
] q = 1, 2, ··· , k −1
(6.66)
然后我们又要构造最高阶的局部截断误差了,另 c
1
= c
2
= ··· = c
p
= 0, 就构造出了 p 阶的多步法,其中
77
数值分析笔记 第六章 常微分方程初值问题数值解法 Tsui Dik Sang
推论 6.4.3 (p 阶多步法).
T
n+k
= c
p+1
h
p+1
y
(p+1)
(x
n
) + O(h
p+2
) (6.67)
其中第一项为局部截断误差,第二项为误差常数
另一方面,相容性的限制为 p ⩾ 1, 就得到了线性多步法相容的充要条件
定理 6.4.4 (线性多步法相容的充要条件). 需要 c
0
= c
1
= 0,即
k−1
X
i=0
α
i
= 1
k−1
X
i=1
iα
i
+
k−1
X
i=0
β
i
= 1
(6.68)
书上在这之后还有对 k=1 的单步法讨论,看看能不能与 Euler 法和梯形法对上,这里就从略了。
6.4.2 Adams 显式与隐式多步法
我们可以将这个方法理解为对于一阶节点是单步的,但是对于一阶导节点是多步的
定义 6.4.3 (Adams 方法).
y
n+1
= y
n
+ h
k−1
X
i=0
β
i
f
n+i
(6.69)
推论 6.4.5 (显式与隐式 Adams 方法). • 显式 Adams 方法:β
k
= 0,
• 隐式 Adams 方法:β
k
̸= 0,
又是要寻找最高阶精度
8
推论 6.4.6 (k=3 时候的 Adams 显式多步法 (3 阶)).
β
0
+ β
1
+ β
2
+ β
3
= 1
2(β
1
+ 2β
2
+ 3β
3
) = 5
3(β
1
+ 4β
2
+ 9β
3
) = 9
4(β
1
+ 8β
2
+ 27β
3
) = 65
(6.70)
不想再说什么,直接打表
8
即满足 c
1
= c
2
= · · · = c
p
= 0
78
数值分析笔记 第六章 常微分方程初值问题数值解法 Tsui Dik Sang
表 6.1: Adams 显式方法(Adams-Bashforth 法)
k p 公式 c
p+1
1 1 y
n+1
= y
n
+ hf
n
1
2
2 2 y
n+1
= y
n
+ h(
3
2
f
n
−
1
2
f
n−1
)
5
12
3 3 y
n+1
= y
n
+ h(
23
12
f
n
−
16
12
f
n−1
+
5
12
f
n−2
)
9
24
4 4 y
n+1
= y
n
+ h(
55
24
f
n
−
59
24
f
n−1
+
37
24
f
n−2
−
9
24
f
n−3
)
251
720
表 6.2: Adams 隐式方法(Adams-Moulton 法)
k p 公式 c
p+1
1 1 y
n+1
= y
n
+ hf
n+1
−
1
2
2 2 y
n+1
= y
n
+ h(
1
2
f
n+1
+
1
2
f
n
) −
1
12
3 3 y
n+1
= y
n
+ h(
5
12
f
n+1
+
8
12
f
n
−
1
12
f
n−1
) −
1
24
4 4 y
n+1
= y
n
+ h(
9
24
f
n+1
+
19
24
f
n
−
5
24
f
n−1
+
1
24
f
n−2
) −
19
720
6.4.3 更多经典线性多步法
6.4.3.1 Milne 法
属于是 k=4 情况下的另外一种最高阶的解, 求解 c
0
= c
1
= c
2
= c
3
= c
4
= 0 的等价条件
β
0
+ β
1
+ β
2
+ β
3
+ = 1
2(β
1
+ 2β
2
+ 3β
3
) = 16
3(β
1
+ 4β
2
+ 9β
3
) = 64
4(β
1
+ 8β
2
+ 27β
3
) = 256
(6.71)
得到
定义 6.4.4 (Milne 法).
y
n+4
= y
n
+
4h
9
(2f
n+3
− f
n+2
+ 2f
n+1
) (6.72)
推论 6.4.7 (Milne 法的局部截断误差).
T
n+4
=
14
45
h
5
y
(5)
(x
n
) + O(h
6
) (6.73)
因此是 4 阶
Milne 法的另外一种推导是通过两端积分下面的式子得到
y(x
n+4
) − y(x
n
) =
ˆ
x
n+4
x
n
f(x, y(x))dx (6.74)
6.4.3.2 Simpson 法
如果积分区间包含的点小一点,编程二步的,就有了 Simpson 法, 即积分
y(x
n+2
) − y(x
n
) =
ˆ
x
n+2
x
n
f(x, y(x))dx (6.75)
79
数值分析笔记 第六章 常微分方程初值问题数值解法 Tsui Dik Sang
定义 6.4.5 (Simpson 法).
y
n+2
= y
n
+
h
3
(f
n
+ 4f
n+1
+ f
n+2
) (6.76)
推论 6.4.8 (Simpson 法的局部截断误差).
T
n+2
= −
h
5
90
y
(5)
(x
n
) + O(h
6
) (6.77)
因此是 4 阶的, 这是在二步方法里面阶数最高的了。
6.4.3.3 Hamming 法
刚刚的 Simpson 是二步方法里面阶数最高的,但是不稳定,因此 Hamming 法追求的事一个稳定性
9
,他是三步的。令 α
1
= 0
去解方程
定义 6.4.6 (Hamming 法).
y
n+3
=
1
8
(9y
n+2
− y
n
) +
3h
8
(f
n+3
+ 2f
n+2
− f
n+1
) (6.78)
推论 6.4.9 (Hamming 法的局部截断误差).
T
n+3
= −
h
5
40
y
(5)
(x
n
) + O(h
6
) (6.79)
6.4.3.4 PECE 预测-校正法
改进的思路和改进 Euler 法类似的,都是为了减少计算量,实际上,改进的 Euler 法就是 PECE 预测-校正法的一种二阶形
式
定义 6.4.7 (PECE 预测-校正法).
预测 P:y
p
n+4
= y
n+3
+
h
24
(55f
n+3
− 59f
n+2
+ 37f
n+1
− 9f
n
)
校正 C:f
P
N
+4
= f (x
n+4
, y
p
n+4
)
校正 C:y
n+4
= y
n+3
+
h
24
(9f
P
n+4
+ 19f
n+3
− 5f
n+2
+ f
n+1
)
校正 C:f
C
N+4
= f (x
n+4
, y
n+4
)
(6.80)
6.4.3.5 PMECME
预测
-
校正法
根据计算我们已经有了表6.1和6.2, 注意到其中的 c
p+1
这一栏,我们完全可以加入在预测之后利用这个修正项去先减小一点
误差, 这也就是这个方法的内容
9
但是至于为什么稳定,书上并没有证明,难道是步数越高越稳定?可能!
80
数值分析笔记 第六章 常微分方程初值问题数值解法 Tsui Dik Sang
定义 6.4.8 (PMECME 预测-校正法).
预测 P:y
p
n+4
= y
n+3
+
h
24
(55f
n+3
− 59f
n+2
+ 37f
n+1
− 9f
n
)
修正 M:y
pm
n+4
= y
p
n+4
+
251
270
(y
c
n+4
− y
p
n+3
)
求值 E:f
pm
n+4
= f (x
n
+4
, y
pm
n+4
)
校正 C:y
c
n+4
= y
n+3
+
h
24
(9f
pm
n+4
+ 19f
n+3
− 5f
n+2
+ f
n+1
)
修正 M:y
n+4
= y
c
n+4
+
19
270
(y
c
n+4
− y
p
n+3
)
求值 E:f
n+4
= f (x
n+4
, y
n+4
)
(6.81)
6.4.3.6 总结
总而言之,还有很多方法,构造完全是见仁见智各取所需,就是精度、计算量以及稳定性之间的权衡了。
6.4.4 收敛性与稳定性
6.4.4.1 相容性
引理 6.4.10 (线性多步法的相容性与收敛性). 收敛则一定相容,相容不一定收敛
定义 6.4.9 (线性多步法的特征多项式). 对于线性多步法
y
n+k
=
k−1
X
i=0
α
i
y
n+i
+ h
k−1
X
i=0
β
i
f
n+i
(6.82)
有
ρ(ξ) = ξ
k
−
k−1
X
i=0
α
i
ξ
i
σ(ξ) =
k
−
1
X
i=0
β
i
ξ
i
(6.83)
称 ρ(ξ) 为第一特征多项式,σ(ξ) 为第二特征多项式
定理 6.4.11 (线性多步法相容的充要条件).
ρ(1) = 0
ρ
′
(1) = σ(1)
(6.84)
后面关于稳定性等的结论相当晕,先从略。
6.4.5 一阶方程组域刚性方程组
6.4.5.1 一阶方程组
即求解微分方程组每一行的元方程为
y
′
i
= f
i
(x, y
1
, y
2
, ··· , y
m
), i = 1, 2, ··· , m (6.85)
81
数值分析笔记 第六章 常微分方程初值问题数值解法 Tsui Dik Sang
令
y =
y
1
y
2
··· y
m
T
f =
f
1
f
2
··· f
m
T
y
0
=
y
0
1
y
0
2
··· y
0
m
T
(6.86)
定义 6.4.10 (一阶方程组).
y
′
= f (x, y)
y(x
0
) = y
0
(6.87)
如果用四阶 Runge-Kutta 法来求解这个方程组,递推公式就是
y
n+1
= y
n
+
h
6
(k
1
+ 2k
2
+ 2k
3
+ k
4
) (6.88)
6.4.5.2 化高阶方程为一阶方程组
先看高阶形式的单个方程
y
(m)
= f (x, y, y
′
, ··· , y
(m−1)
)
y(x
0
) = y
0
y
′
(x
0
) = y
1
.
.
.
y
(m−1)
(x
0
) = y
m−1
(6.89)
则可以设置高阶为多个未知量,即
y
1
= y
y
2
= y
′
.
.
.
y
m
= y
(m−1)
(6.90)
则高阶方程就转化成了一阶方程组
6.4.5.3 刚性方程组
定义 6.4.11 (刚性问题). 解的分量在数量级上相差很大
下面不加证明的给一个结论
定理 6.4.12 (一般线性系统的通解). 对于线性系统
dy
dt
= Ay(t) + g(t) (6.91)
其中
y, g ∈ R
N
A ∈ R
N×N
(6.92)
82
数值分析笔记 第六章 常微分方程初值问题数值解法 Tsui Dik Sang
若 A 的特征值 λ
j
= α
j
+ iβ
j
相应的特征向量为 ψ
j
, 则该线性系统的通解为
y(t) =
N
X
j=1
c
j
e
α
j
t
φ
j
+ ψ(t) (6.93)
其中 ψ(t) 为特解
推论 6.4.13 (稳态解). 如果 Re(λ
j
) < 0, 带入易得
lim
t→∞
y(t) = ψ(t) (6.94)
此时 ψ(t) 为稳态解
所以特征值实部小于零是稳态的前提条件,之后可以再定义反应刚性程度的系数
定义 6.4.12 (刚性系数). 若线性系统方程 A 的特征值均满足实部小于零, 且
s =
max
1⩽j⩽N
|Re(λ
j
)|
min
1⩽j⩽N
|Re(λ
j
)|
>> 1 (6.95)
则称这个方程为刚性方程,s 称为刚性系数
我们在这门课中提及这个概念也只是一个定性的了解的,Gear 方法等可以在求解刚性方程组的时候使用,不再介绍
83
数值分析笔记 第六章 常微分方程初值问题数值解法 Tsui Dik Sang
84
第七章 笔记初稿完成
至此,初稿笔记就将考试内容的所有内容就整理完了。
Tsui Dik Sang
—2025.5.5
结合复习过程又补充了一些遗漏点
Tsui Dik Sang
—2025.5.15
85
数值分析笔记 第七章 笔记初稿完成 Tsui Dik Sang
86
参考文献
[1] 聂玉峰, 王振海. 数值方法简明教程(第二版). 高等教育出版社, 2020. ISBN 9787040552102.
[2] 李庆扬. 数值分析(第 5 版). 高等教育出版社, 2008. ISBN 9787302185659.
87