当我学数值分析时,我学到了什么?
learning why, thinking what, then forgetting how.
随着时间的流逝,知识总会被遗忘和被沉淀,我们无法选择去遗忘那一部分,但是我们可以选择去沉淀那一部分。
本人配套开发的,基于 Go 的线性代数库 Singular 已开源,实现了基础的数值分析算法。
教材为:《数值分析 第4版(颜庆津 编著)》
第一章 绪论
研究对象
科学计算:通过建立数学模型把科学技术问题转化为数学问题,然后对数学问题进行离散化,将其转化为数值问题,最后使用数值计算方法计算出数值问题的解。
误差和算法
误差的来源与分类:
- 模型误差
- 观测误差
- 截断误差
- 舍入误差
在数值分析中,主要研究截断误差和舍入误差对计算结果的影响,而一般不考虑模型误差和观测误差。
函数求值的误差估计:由参数值的误差可以估计出函数值的误差
设计算法的原则:
- 要有数值稳定性,即能控制舍入误差的传播。
- 即在迭代中舍入误差的系数应小于一,事实上,如果大于一,稍作变换处理即可使系数小于一
- 防止小数加到大数
- 防止相近的数相减
- 防止除数绝对值远小于被除数
向量范数与矩阵范数
提出了向量范数的概念,需要满足范数的三个性质:
- 正定性;
- 齐次性
- 三角不等式
提出了矩阵的向量范数的概念,同样需要满足范数的三个性质。如果同时满足相容性,则成为矩阵的矩阵范数。
向量范数和矩阵范数也可以满足相容性。对一个矩阵范数,必存在与之相容的向量范数。对一个向量范数,可以构造与之相容的矩阵范数,此时称为矩阵的算子范数。
第二章 线性方程组的解法
问题描述:求解线性方程组:AX = b
Gauss 消去法
- 顺序 Gauss 消去法
- 消元;回代
- 列主元素 Gauss 消去法
- 消元(行交换);回代
列主元素 Gauss 消去法的精度更高:如果顺序 Gauss 消去法的主元素绝对值较小,会使得行乘数的绝对值较大,导致小数加到较大数,并积累舍入误差。
直接三角分解法
矩阵的 A = LU 分解不一定存在,有可能为 QA = LU 分解。
- Doolittle 分解:单位下三角矩阵和上三角矩阵
- Crout 分解:下三角矩阵和单位上三角矩阵
求 LU 分解有直接解法。
选主元的三角分解法:QAX = LUX = Qb,回代:
- UX = y
- Ly = Qb
矩阵的条件数与病态线性方程组
研究一下:当 A 和 b 有了微小变化之后,如何影响解向量 x 的变化。可以由此判断 A 和 b 的误差对 x 的误差的影响。
提出了条件数的概念:
- c o n d ( A ) = ∣ ∣ A ∣ ∣ ⋅ ∣ ∣ A − 1 ∣ ∣ cond(A) = ||A||·||A^{-1}|| cond(A)=∣∣A∣∣⋅∣∣A−1∣∣
- c o n d ( A ) 2 = ∣ λ m a x λ m i n ∣ cond(A)_2 = |\frac{λ_max}{λ_min}| cond(A)2=∣λminλmax∣
条件数较大时,即为病态方程组,即当 A 和 b 有微小变化时,原方程组的解会有很大的相对误差。
优化病态方程组的求解:
- 采用高精度的算数运算。例如双精度。
- 平衡方法。使用列平衡或行平衡来降低条件数:DAx = Db。
- 残差校正法。逐步迭代,累加解向量的修正量。
迭代法
A x = b Ax = b Ax=b:
- A = N − P A = N - P A=N−P
- N x = P x + b Nx = Px + b Nx=Px+b
- x = G x + d , G = N − 1 P , d = N − 1 b x = Gx + d, G = N^{-1}P, d = N^{-1}b x=Gx+d,G=N−1P,d=N−1b
分解系数矩阵,构造不动点方程,迭代求解解向量。
迭代收敛的充分必要条件是 G 的谱半径小于一。
分解系数矩阵为上三角矩阵,对角矩阵,下三角矩阵,A = D + L + U:
- Jacobi 迭代法: N = D , P = − ( L + U ) N = D, P = -(L + U) N=D,P=−(L+U)
- Gauss-Seidel 迭代法: N = D + L , P = − U N = D + L, P = -U N=D+L,P=−U
- 逐次超松弛迭代法(SOR): N = 1 w D + L , P = − [ ( 1 − 1 w ) D + U ] N = \frac{1}{w}D + L, P = -[(1-\frac{1}{w})D + U] N=w1D+L,P=−[(1−w1)D+U]
Jacobi 迭代法和 Gauss-Seidel 迭代法的比较:其实高斯-赛德尔迭代和雅克比迭代相比,只是用到了新计算的值进行迭代(从分量写法可以看出来),但是两者的迭代速度并没有好坏之分。Jacobi 迭代法的迭代矩阵保持稀疏性,即 A 稀疏 -> G 稀疏;而 Gauss-Seidel 迭代法的迭代矩阵不保持稀疏性。保持稀疏性可以使得计算量减少。而且 Gauss-Seidel 迭代法由于使用了新计算的值,所以无法并行求解。
Gauss-Seidel 方法是 SOR 方法的特殊情况(w=1)。SOR 方法中,调整 w,使得 x 的迭代序列变得稀疏(松弛),速度收敛就可以加快。
第三章 矩阵特征值与特征向量的计算
幂法和反幂法
- 幂法:迭代求解按模最大特征值及其特征向量
- 幂法实质是构造了不动点方程迭代求解
- 反幂法:迭代求解按模最小特征值及其特征向量
- 反幂法实质是求解方程组求解逆矩阵的按模最大特征值
- 原点位移:由反幂法求解距离位移按模最近的特征值及其特征向量
- 利用了矩阵的原点平移变换,特征值变换相同
Jacobi 方法
根据正规分解,迭代求实对称矩阵(正规矩阵的特殊形式)的全部特征值和特征向量。
QR 方法
求 QR 分解有逐步 Householder 解法。
QR 方法可以求解全部特征值,其实质是将矩阵由 Q 相似变换到拟上三角阵。相似变换保特征值不变,且上三角阵的特征值在其对角线上,当拟上三角阵逼近上三角阵时,其对角线即为特征值的逼近。
- QR 方法迭代求解。
- 带双步位移的 QR 方法迭代求解。
- 可以减少迭代次数,但是每次迭代的运算增加(2 次 QR);还可以进行修正,使得计算量减少(1 次 QR):
第四章 非线性方程与非线性方程组的迭代解法
问题描述:求解非线性方程:f(x) = b
基础理论
- 大范围收敛定理:在区间内任取初值,迭代都收敛
- 局部收敛定理:在接近根的邻域内取初值,迭代收敛
收敛速度的衡量: l i m k − > + o o ∣ e k + 1 ∣ ∣ e k ∣ r ≤ c lim_{k -> +oo} \frac{|e_{k+1}|}{|e_k|^r} ≤ c limk−>+oo∣ek∣r∣ek+1∣≤c,则收敛速度为 r 阶。r 越大收敛越快,表示迭代一步,误差变小了很多。
- r = 1:线性收敛
- r = 2:平方收敛
- r > 1:超线性收敛
非线性方程的迭代解法
- 对分法
- 只能求解单根,奇数重根。
- 简单迭代法(不动点迭代法)
- 构造不动点方程。
- 若不动点方程的前 m 阶导数为 0,则收敛阶为 m+1。
- Steffensen 迭代法
- 在简单迭代法的不动点方程导数不为零时,为线性收敛。
- 使用 Steffensen 迭代法对不动点方程进行构造后,可以使达到平方收敛。
- 但是如果原不动点方程的已有大于2阶的收敛速度,则使用 Steffensen 构造对提高收敛速度的作用不大。
- Newton 法(切线法)
- 使用一阶导数,类似梯度下降法。
- 至少是平方收敛。
- 求 m 重根的 Newton 法
- 稍微修改下 Newton 法的迭代函数即可。
- 只有线性收敛速度。
- 当重根数 m 未知时,可以进行函数转化为求新函数的单零点。缺点是需要计算二阶导数。
- 割线法
- 使用增量比(割线)替换 Newton 法的导数(切线)。
- 不需要求解导数,但是收敛速度低于 Newton 法,高于线性收敛,收敛阶至少为 1.618。
- 单点割线法
- 在割线法中固定开始点。
- 属于简单迭代法,线性收敛速度。
非线性方程组的迭代解法
- 简单迭代法(不动点迭代法)
- 线性收敛速度。
- Newton 法
- 至少平方收敛速度,需要求解 Jacobi 矩阵。
- 一般对初始向量要求接近解向量,即局部收敛。
- 离散 Newton 法
- 使用差商代替求导,不需要求解 Jacobi 矩阵。
- 一般收敛阶低于 Newton 法,如果合理选择步长,即为 Newton-Steffensen 法,也可以达到至少平方收敛速度。
第五章 插值与逼近
插值
问题描述:求解插值函数: f ( x ) = p n ( x ) = ∑ k = 0 n c k φ k ( x ) f(x) = p_n(x) = \sum^n_{k=0} c_k φ_k(x) f(x)=pn(x)=∑k=0nckφk(x)
首先提出了函数组线性无关的概念,即函数值构成的向量线性无关。
满足条件的 n 次插值多项式是存在且唯一的。
代数插值
只要求插值函数的函数值。
- Lagrange 插值法:使用 Lagrange 插值多项式
- Newton 插值法:使用 Newton 插值多项式
n+1 个插值节点的 n 次插值多项式是存在且唯一的。
插值多项式的误差主要在于截断误差。其影响因素在于插值多项式的次数和插值节点的选择。
理论和实践中:高次插值不可取,通常使用分段低次插值。不计算整个插值函数,而是针对每个x选择最靠近x的若干个(2~4)节点作为插值节点,计算分段插值函数。不适用于光滑性要求高的外形设计。
Hermite 插值
同时要求插值函数的函数值和导数值。
n+1个函数值,m+1个导数值,其n+m+1次插值多项式是存在且唯一的。
- 先求得函数值的插值多项式(Lagrange 或 Newton 插值法)。
- 加上函数值为零的尾项使其满足导数值条件。
样条插值
使用分段低次多项式插值,且能满足光滑性要求。
要求各个节点处的被插值函数值和两个边界节点的导数值。
边界条件:
- 边界节点处的二阶导数相等
- 边界节点处的一阶导数相等
- 周期函数:两个边界节点的一阶导数相等
在实际应用中,不仅用三次样条函数 s ( x ) s(x) s(x)计算被插值函数 f ( x ) f(x) f(x)的近似值,而且常用 s ′ ( x ) s^{'}(x) s′(x)近似计算 f ′ ( x ) f^{'}(x) f′(x)。
三角插值
三角插值等于傅里叶变换。
快速傅里叶变换: N 2 − > N l o g N N^2 -> NlogN N2−>NlogN
正交多项式
在[a, b]上带权ρ(x)的正交函数系一定是在[a, b]上线性无关的函数系,而不论ρ(x)是什么。
可以使用 Gram-Schmidt 正交化方法构造正交多项式系。
常用的正交多项式:
- Legendre 多项式
- Chebyshev 多项式
- Laguerre 多项式
- Hermite 多项式
最佳平方逼近
给定子空间,则对于f(x)的最佳平方逼近元素是唯一的。可用法方程或正规方程求最佳逼近平方元素,为了避免系数矩阵的是病态矩阵,应采用正交基底。选择不同的基函数,可以求得不同的最佳平方逼近。
拟合不要求过点,而插值要求过点。
第六章 数值积分
计算定积分只要求出被积函数的原函数,再利用牛顿-莱布尼兹公式即可。但是有很多被积函数是很难求出它的原函数的,甚至不能求出有限形式的原函数。即无法采用解析方法计算定积分,本章介绍数值积分法。
∫ a b f ( x ) d x ≈ ∑ k = 0 n λ k f ( x k ) \int^b_a f(x) dx ≈ \sum^n_{k=0} λ_k f(x_k) ∫abf(x)dx≈k=0∑nλkf(xk)
代数精确度:对求积公式,当f(x)为任何次数不高于m的多项式时都成为等式,而当为次数m+1次多项式时不能成为等式,则称它具有m次代数精确度。
求解:forgetting how
第七章 常微分方程初值问题的数值解法
求 解 : y ( t ) , { y ′ = f ( t , y ) y ( t 0 ) = y 0 求解:y(t), \left\lbrace \begin{array}{c} y^{'} = f(t, y)\\ y(t_0) = y_0 \end{array} \right. 求解:y(t),{y′=f(t,y)y(t0)=y0
欧拉方法分为显式积分和隐式积分两种形式,其中显式积分条件稳定,隐式积分无条件稳定。
具体实际问题应用可以参考:显式积分,隐式积分和弹簧质点系统
求解:forgetting how
第八章 偏微分方程的差分解法
求解泊松方程: ∂ 2 u ∂ x 2 + ∂ 2 v ∂ y 2 = f ( x , y ) \frac {\partial ^2 u} {\partial x^2} + \frac {\partial ^2 v} {\partial y^2} = f(x,y) ∂x2∂2u+∂y2∂2v=f(x,y)
具体实际问题应用可以参考:从泊松方程的解法,聊到泊松图像融合
求解:forgetting how