《数值分析》笔记

使用教材:《数值分析》李庆扬,王能超,易大义 编。
大致梳理整本书的内容和重点。

第1章 数值分析与科学计算引论

本章是介绍性质的,主要说了数值分析(乃至计算数学)在学什么,解决什么问题。
内容上,需要掌握误差的几个种类,绝对/相对误差(限)的计算,有效数字。
1.3节“误差定性分析与避免误差危害”:了解数值稳定、病态问题、条件数的概念。尤其是条件数的意思,在后续数值线性代数中也会用到。避免误差危害的方法:避免两相近数相减,避免用绝对值小的数做除数,注意运算次序和减少运算次数。
了解一些古典的数值计算思想:秦九韶算法,迭代法,以直代曲,加权平均的松弛技术(感觉类似外推法?)
本章内容浅显零散,了解即可。

第2章 插值法

注意区分插值拟合:插值函数一定确保已有的数据点完全拟合,一般维数(比如多项式次数)较高;拟合则不需要,一般维数较低。

2.1 引言

差值问题的提出:略
多项式拟合:对于n+1个数据点(x,y),用多项式拟合,则使用n次多项式刚好保证存在唯一。(代数方法证明)
虽然结果唯一,但插值多项式系数的求法有许多种,这里介绍拉格朗日插值牛顿插值

2.2 拉格朗日插值

思路:基函数l_k(x)(线性插值基函数):满足l_j(x_j)=1,l_j(x_i)=0,i\neq j的n+1次多项式(显然唯一)。插值多项式L(x)=\sum y_kl_k(x)为基函数的线性组合,系数由数据点确定。
基函数方法是将插值问题划归于特定条件下容易实现的插值问题,本质上是广义的坐标系方法。
插值余项:若f^{(n+1)}(x)存在,则R_n(x)=f(x)-L(x)=\frac{f^{(x+1)(\xi)}}{(n+1)!}\prod_{i=0}^n (x-x_i)
使用中值定理证明。(中值定理和泰勒展开在数值分析中用处很多)。由插值余项的公式得,要保证多项式插值与原函数较为接近,最好原函数的n+1次导数有界。
拉格朗日插值法思路清晰,公式漂亮,但若数据点变化需要重新计算。下面的牛顿插值多项式解决了这个问题。

2.3 均差与牛顿插值多项式

引入“均差”(差商)概念,类比导数。
牛顿法思路:利用均差,依次添加数据点,每添加一个数据点插值多项式P(x)会增加一项。
对数据点的顺序没有要求。
牛顿插值法在添加数据点时计算量较小,但公式没那么漂亮。

2.4 埃尔米特插值

拉格朗日插值和牛顿插值仅保证数据点f(x_i)=y_i,没有导数信息。埃尔米特插值用于一些f^{(k)}(x_i)=y_i^k的需求。
计算方法类似牛顿插值,必要时可用待定系数法。
书中仅讨论了两个典型的埃尔米特插值的计算。
注意两点三次埃尔米特插值三次样条插值的区别:两者都是用三次多项式对两点做插值,但前者是为了保证f(x_1),f(x_2),f'(x_1),f'(x_2)均与数据相同,后者仅保证f(x_1),f(x_2)与数据相同,导数的条件是为了保证f的光滑性。

2.5 分段低次插值

了解高次插值的病态性质:龙格现象,高次插值没有收敛性。
通过分段,降低次数,减少病态现象。
常用的:分段线性插值,分段三次埃尔米特插值,三次样条插值。
分段插值具有一致收敛性。

2.6 三次样条插值

保证数据点的同时插值函数具有二阶光滑性,在很多实际问题中有应用。
思路:对于n+1个数据点,分成n段三次多项式,共需4n个参数。限制条件:每个数据点吻合,2n个方程;为保证光滑性,每个连接点两边的一阶导数、二阶导数相同,2(n-1)个方程;还差2个方程,一般为两个端点的条件(边界条件),可根据实际问题要求给定。
到此为止已经列出了一个线性方程组,解法看书即可,最终化成三对角矩阵的线性方程组,使用追赶法求解。
性质:插值函数及其一阶导数、二阶导数均一致收敛与数据的原函数。

第3章 函数逼近与快速傅里叶变换

3.1 函数逼近的基本概念

函数逼近问题就是:在区间[a,b]上用简单函数逼近已知复杂函数的问题,简单函数一般是n次多项式、有理函数、分段低次多项式等。(可以理解为,从任意函数构成的空间到多项式函数空间的投影?)
维尔斯特拉斯定理:闭区间上连续函数,均可用有限维的多项式任意逼近(一致范数)。证明方法用到伯恩斯坦多项式,这是一个理论上对[0,1]上连续函数f任意逼近的多项式,但实际上收敛速度慢。

介绍范数与赋范线性空间,用来描述逼近的程度。引入内积、正交、施瓦茨不等式、权函数、函数内积、函数范数等概念。

最佳逼近:对于某个函数空间H,某个范数,找出H中与目标函数f的差的范数最小的函数P*,为最佳逼近函数。
||f(x)-P*(x)||=min_{P\in H_n}||f(x)-P(x)||

若取P为多项式,范数为一致范数,则为最佳一致逼近多项式;若范数取二范数,则为最佳平方逼近多项式;若f(x)为列表函数(离散情况),范数取二范数,则为最小二乘拟合

3.2 正交多项式

正交多项式是函数逼近的重要工具,在数值积分中也有重要应用。
正交多项式可以理解为多项式内积空间中的正交基?

概念:带权\rho (x)正交,正交函数族,标准正交函数族,施密特正交化,正交多项式的递推关系式,正交多项式零点的定理。

勒让德多项式P_n(x):区间为[-1,1],权函数\rho \equiv1时,由\{1,x, \cdots ,x^n,\cdots \}\正交化得到的多项式。
性质:正交性,奇偶性,递推关系,n个不同零点。

切比雪夫多项式T_n(x):区间为[-1,1],权函数\rho = \frac{1}{\sqrt{1-x^2}}时,由\{1,x, \cdots ,x^n,\cdots \}\正交化得到的多项式。另一个比较直观的表示:若令x=cos\theta,则T_n(x)=cos(n\theta),0 \le \theta\ \le \pi
性质:递推关系,(带权)正交性,奇偶性,n个零点,首项系数,最佳一致逼近多项式,切比雪夫多项式零点插值(通过改变插值点的选位,从等距点改为切比雪夫多项式零点,可避免龙格现象)

一些其他有用的正交多项式:第二类切比雪夫多项式,盖拉尔多项式,埃尔米特多项式。

正交函数族的优势在下面一节中就有体现。

3.3 最佳平方逼近

数学表达式:
||f(x)-S*(x)||_2^2=min_{S(x)\in \varphi} ||f(x)-S(x)||_2^2
连续形式:min_{S(x)\in \varphi} \int_a^b \rho(x)[f(x)-S(x)]^2dx
对多项式空间,等价于求多元函数的最小值:I(a_0,a_1,\cdots,a_n) = \int_a^b \rho(x) \left [ \sum_{j=0}^n a_j \varphi_j(x) - f(x) \right ]^2 dx
这是一个分析性质比较好的函数,求最小值只要求某个点对所有变量的偏导数等于0,刚好是一个线性方程组,称为法方程。这样看起来是找到了解析解,只要用数值方法求解线性方程组就可以,但直接取\varphi_k(x)=x^k时,系数矩阵为希尔伯特矩阵,是一个条件数大的矩阵,对它解线性方程组得到的误差较大。

为了解决上面的问题,可以用正交函数族作最佳平方逼近,参考傅里叶级数,由正交性质,系数矩阵是对角阵,因此在计算傅里叶级数的系数时不需要解线性方程组。在上面的例子中将x^k换成勒让德多项式P_k(x)就能解决问题。这就体现出正交函数族对于数值逼近的优势。

切比雪夫级数:将函数f按切比雪夫多项式展开,得到切比雪夫级数,可作为f[-1,1]上的近似最佳一致逼近多项式。

3.4 曲线拟合的最小二乘法

类似于最佳平方逼近的离散形式,在科学实验中对实验数据的处理(曲线拟合)经常用到。
数学表达式:min_{S(x)\in \varphi } \sum_{i=0}^m [S(x_i) - y_i]^2
方法与最佳平方逼近类似,法方程,解线性方程组,用正交多项式解决病态系数矩阵问题。略。

3.5 有理逼近

前面的多项式逼近对于有界连续函数逼近效果较好,现在对于在某点附近无界的情况,用有理函数逼近效果较好。
有理函数:R_{nm}(x)=\frac{P_n(x)}{Q_m(x)} = \frac{\sum_{k=0}^n a_kx^k}{\sum_{k=0}^m b_kx^k},即两个多项式的比。这样可以得到比如\frac{1}{x}这样的可以趋于无穷的函数逼近。
计算方法是利用连分式。

帕德逼近:如果有理函数R_{nm}(x)逼近f(x)R_{nm}^{(k)}(0)=f^{(k)}(0),k=0,1,\cdots,(n+m),则称R_{nm}(x)f(x)x=0处的(n,m)阶帕德逼近。这个逼近的好处个人没用过不太清楚,略。

3.6 三角多项式逼近与快速傅里叶变换

在模型数据具有周期性时,用三角函数特别是正弦函数和余弦函数作为基函数是合适的。
最佳平方三角逼近:傅里叶级数
快速傅里叶变换:略。

评注

正交多项式中的勒让德多项式和切比雪夫多项式是两个十分重要且经常使用的正交多项式,应引起高度关注。
当模型为多项式时法方程是病态的,为此推荐用正交化方法避免解法方程。
如果数据是周期的,使用三角最小二乘或三角插值是合适的,计算用快速傅里叶变换(FFT),它是节省计算量的一个范例。

第4章 数值积分与数值微分

4.1 数值积分概论

一般的数值积分:\int_a^b f(x) dx \approx \sum_{k=0}^n A_k f(x_k)
其中x_k为求积节点,A_k为求积系数,即节点x_k的权重。
代数精度:描述求积公式的误差的阶数,计算方法为计算f(x) = x^k的情况的误差。
插值型求积公式:先求f(x)的插值多项式L_n(x),再用I_n := \int_a^b L_n(x) dx作为I := \int_a^b f(x) dx的近似值,其误差可用插值多项式的误差计算:R[f] = \int_a^b [ f(x) - L_n(x) ] dx = \int_a^b R_n(x) dx
求积公式是插值型的,当且仅当其至少有n次代数精度。

收敛性:只要插值点充分密集,数值积分就会充分接近真实值。
稳定性:插值点数据(x_k,f(x_k))可能存在误差,稳定性是指,只要误差充分小,数值积分就会充分接近真实值。反之,称为对误差敏感。

4.2 牛顿-柯特斯公式

先默认选取等距节点,因为简单,此时的公式I_n = (b-a) \sum_{k=0}^n C_k^{(n)} f(x_k)称为牛顿-柯特斯公式,C_k^{(n)}称为柯特斯系数。

对于[a,b],选取a、b两点插值(n=1),称为梯形公式;选取a、b、(a+b)/2三点(n=3),称为辛普森公式,5个点(n=4)称为柯特斯公式。9点以上的牛顿柯特斯公式稳定性差,不用。

当n为偶数,牛顿柯特斯公式至少具有n+1阶代数精度。

4.3 复合求积公式

类似分段低阶插值,先分段再用低阶求积公式,称为复合求积公式。常用的,复合梯形公式(折线段),复合辛普森公式。

4.4 龙贝格求积公式

思路:实际计算时若精度不够可将步长逐次分半,直到精度足够为止。

理查森外推法:只要真值与近似值的误差能表示成h的幂级数,就能使用外推法提高精度。
例如,T(h) = I + a_1 h^2 + a_2 h^4 + \cdots,则可增加一些点将h减半得到T(h/2) = I + \frac{1}{4}a_1 h^2 + \frac{1}{16}a_2 h^4 + \cdots,通过线性组合得到S(h) = \frac{4}{3} T(h/2) - \frac{1}{3} T(h) = I + b_1 h^4 + b_2 h^6 + \cdots,误差阶数就从2阶提高到了4阶!

龙贝格算法:从复合梯形公式开始,不断使用外推法,(计算有递推公式),直到精度足够为止。

4.5 自适应积分方法

回顾以上符合积分公式,对整个积分区域使用同样的阶数,适用于被积函数变化不大的积分。
对于在不同区域变化程度不同的函数,我们可以对不同区间使用不同的插值点数量(根据误差自动判断,龙贝格算法),平坦的部分步长大,变化剧烈的部分步长小,既减少运算量又减少误差。具体过程略。

4.6 高斯求积公式

以上我们都默认等距节点,现在我们改变节点的位置,希望在不改变节点数量的情况下得到更高阶的误差,这就是高斯求积公式。
如果求积公式(n+1个节点)具有2n+1次代数精度(n+1个节点的理论最高精度),则称其节点为高斯点,公式为高斯型求积公式。

如何取高斯点?
由定理证明。。。高斯点是正交多项式的零点。
性质:高斯求积公式收敛,稳定。

高斯-勒让德求积公式:使用勒让德多项式作为高斯求积公式需要的正交多项式。
对于区间不是[-1,1]的情况,做变换x=\frac{b-a}{2}t + \frac{a+b}{2}即可。

高斯-切比雪夫求积公式:与上面类似,应对不同情况。
还有其他正交多项式。

注:高斯求积公式不一定要用到区间的端点,对一些在端点出趋于无穷的函数积分(反常积分)近似较好。

4.7 多重积分

写成累次积分,化成一重积分的情况。略。

4.8 数值微分

思路1:用差商近似导数,比如中点公式f'(a) \approx \frac{f(a+h) - f(a-h)}{2h},误差用泰勒展开计算。

思路2:插值型求导公式:先求f的插值函数P_n(x),再用P'_n(x)作为导数的近似值。

思路3:三次样条求导。利用三次样条函数导数值也收敛的性质,可作为数值微分的方法。

思路4:数值微分的外推算法,对思路1的改进,略。

第5章 解线性方程组的直接方法

5.1 引言与预备知识

线性方程组的系数矩阵可大致分为:低阶稠密矩阵,大型稀疏矩阵(零元素多)。
方法大致分两种:直接法(对低阶稠密矩阵有效),迭代法(对大型稀疏矩阵有效)。(本章只讲前者)
特殊形式的矩阵(三角,三对角等)有特殊方法。

概念:矩阵特征值,谱半径,特殊矩阵,若当型矩阵

5.2 高斯消去法

分为前代法和回代法两步,分别求解三角矩阵方程,LU分解(A=LU)。具体可参考《数值线性代数》。

约化的主元素都不为零(普通高斯消去法能进行下去)的充要条件是,矩阵A的顺序主子式都不为零。

对于约会的主元有零的情况:列主元消去法,每次循环时增加一步行变换,可避免主元是零或接近零,减少误差。只要矩阵可逆就能计算下去,实用。列主元的LU分解:PA=LU

5.3 矩阵三角分解法

普通的LU分解和列主元LU分解:上面提到了,具体步骤略。

平方根法:对于对称正定矩阵的专用方法。A=LDL^T=L_1L_1^T楚列斯基分解)。具体步骤略。

追赶法:对于三对角矩阵的专用方法。具体步骤略。

特殊矩阵的专用方法比平凡的方法效果更好(时间、空间复杂度低),但适用范围小。

5.4 向量和矩阵范数

度量,为矩阵/向量的误差分析做铺垫。

矩阵的算子范数/诱导范数:由向量范数扩展到矩阵,具有相容性。也有不是算子范数的矩阵范数,比如F范数(弗罗贝尼乌斯范数)。

有一些定理需要了解。

5.5 误差分析

病态方程组、病态矩阵:对矩阵A或常数项b的变化敏感。条件数刻画了病态程度。希尔伯特矩阵是病态矩阵。若条件数计算不方便,也有几条判断病态的简单方法。

面对病态问题,首先是列主元方法,若解决不了还可以采用高精度算术运算或采用预处理方法

第6章 解线性方程组的迭代法

6.1 迭代法的基本概念

迭代法的格式:x^{(k+1)} = Bx^{(k)} + f
先(任意)取一个初值x^{(0)},再有Ax=b的A和b确定迭代格式,不断迭代,直到达到一定次数或误差达到要求为止。
如果\lim _{k \rightarrow \infty } x^{(k)}存在,则称迭代法收敛。

定义向量序列的极限和矩阵序列的极限。
\lim _{k \rightarrow \infty }B^k=0等价于\rho(B)<1等价于B的某个范数||B||<1

对于任意初值x^{(0)},迭代法x^{(k+1)} = Bx^{(k)} + f收敛的充要条件是\lim _{k \rightarrow \infty }B^k=0

迭代法收敛速度的度量:平均收敛速度,渐进收敛速度。

6.2 雅克比迭代法与高斯-赛德尔迭代法

对于Ax=b,令A=D-L-U,D为对角矩阵,L为下三交矩阵,U为上三角矩阵。
则有Dx = (L+U)x+b,即x = D^{-1}(L+U)x+D^{-1}b,此为雅克比迭代法的迭代格式。
用方程组表示,就是给定x^{(k)},对每个元素x^{(k+1)}(i)迭代时都使用x^{(k)}的元素,全部算完再更新为x^{(k+1)}

类似的,有(D-L)x = Ux+b,即x = (D-L)^{-1}Ux+(D-L)^{-1}b,此为高斯-赛德尔迭代法的迭代格式。
用方程组表示,就是计算每个元素x^{(k+1)}(i)尽可能使用最新的结果,即x^{k+1}(j)x^{k}(l)\forall j<i, \forall l >i

收敛性:(可用迭代法收敛条件判断)
概念:对角占优矩阵,严格对角占优矩阵,可约矩阵,不可约矩阵。
对角占优定理:如果A为严格对角矩阵或不可约对角占优矩阵,则A为非奇异矩阵。
如果A为严格对角矩阵或不可约对角占优矩阵,雅克比迭代法和高斯赛德尔迭代法均收敛。

6.3 超松弛迭代法

高斯赛德尔迭代法的改进,SOR超松弛迭代法。略。

6.4 共轭梯度法

很重要!!需要仔细看书,在此就不细说了,只写需要掌握的部分。

对于对称正定矩阵A,考虑二次函数\varphi:R^n \rightarrow R
\varphi(x) = \frac{1}{2}(Ax,x) - (b,x) = \frac{1}{2}\sum_{i=1}^n \sum_{j=1}^n a_{ij}x_ix_j - \sum_{j=1}^n b_j x_j
\varphi(x)是可导的凸函数,其去最值的时候等价于梯度等于零,即\bigtriangledown \varphi(x) = Ax-b=0
因此求解线性方程组的问题等价于求二次函数最小值的问题,求解方法是构造一个向量序列\{ x^{(k)} \}\使二次函数趋向于最小值。

最速下降法:每次迭代方向为\varphi(x)x^{(k)}处的梯度方向,确定方向后取**该方向上\varphi(x)的最小值点为迭代结果x^{(k+1)}。即以梯度方向为迭代方向,再用线搜索确定步长。
计算:梯度:先求出梯度函数即可。线搜索:令导数等于零可求得步长。
性质:两个相邻的搜索方向正交,因此快接近最优值时会出现锯齿状的路线,速度变慢。收敛速度与矩阵A的最大最小特征值的差有关,特征值越接近收敛速度越快。
速度并不快的方法,但理论简单,后续的方法都是在其上的改进。

共轭梯度法(CG方法):求解大型稀疏对称正定方程组十分有效的方法。
思路:搜索方向从单纯的梯度方向改进为,x^{(k)}处梯度方向和上一步的收缩方向组成的平面,这样需要两个参数来确定步长。
性质:误差r^{(k)}:=Ax^{(k)}-b序列是正交向量组,搜索方向p^{(k)}序列是A-共轭向量组。
r^{(k)}正交性,对n维线性方程组,理论上最多n步便可得到精确解,可以说CG方法是直接法;但由于舍入误差,很难保证r^{(k)}的正交性,一般还是以误差到达精度要求时终止,也可以说CG方法是迭代法。

共轭梯度法的细节可参考《数值线性代数》,它还有一些其他的性质。
目前讨论的共轭梯度法是关于对称正定矩阵A的,对于一般矩阵也可以稍加修改后使用,具体参考《最优化理论与方法》。

第7章 非线性方程与方程组的数值解法

7.1 方程求根与二分法

二分法:最朴素的方法,在函数连续且端点函数值异号的情况下,不断压缩区间长度直到达到精度要求为止。

7.2 不动点迭代法及其收敛性。

x = \varphi(x),则称x\varphi(x)的不动点。若迭代法收敛,则最终收敛到其不动点。

不动点的存在性与迭代法收敛性:不动点存在唯一性定理
\varphi(x)[a,b]连续,且a \le \varphi(x) \le b,且存在正常数L<1使\forall x,y \in [a,b]. |\varphi(x) - \varphi(y)| \le L |x-y|,则\varphi(x)[a,b]上存在唯一的不动点x^*

局部收敛性:若varphi(x)有不动点x^*,如果存在不动点附近的邻域,使得当初值在此邻域中时迭代法一定收敛到不动点,则称迭代法局部收敛
局部收敛的充分条件:\varphi'(x)在不动点附近连续且绝对值小于1,则局部收敛。

收敛速度:迭代误差序列(一个无穷小量)的阶数,p阶收敛。

7.3 迭代收敛的加速方法

已有迭代格式\varphi,通过一些变形得到新的迭代格式\phi[\varphi],使新的迭代格式收敛速度(阶数)加快。
埃特金加速收敛方法,斯特芬森迭代法。

7.4 牛顿法

思路:利用导数(梯度)信息,对原函数线性化(泰勒展开)来近似,再求零点。
x_{k+1} = x_k - \frac{f(x_k}{f'(x_k)}
有几何解释:切线法。
至少是二阶收敛,优点是收敛速度快,缺点是每一步计算量稍大(要计算导数值),而且对初值有要求(不是全局收敛,需要初值在不动点附近)
可以扩展到高维情况,为了克服缺点有一些变种(简化牛顿法,牛顿下山法),重根情形可通过变形加速。

7.5 弦截法与抛物线法

以上是用前一个点迭代后一个点,称为单点迭代法。
弦截法与抛物线法使用前几个点计算后一个点,称为多点迭代法,从而可以不使用导数信息且保持高阶收敛速度。

7.6 求根问题的敏感性与多项式的零点

知道敏感性的意思即可。

求多项式的全部零点:以上方法只能求一个零点。可以每求一个零点,利用因式分解给原方程降次。推荐使用抛物线法,可求复根。

7.7 非线性方程组的数值解法

在一维的基础上向高维推广。使用向量范数代替误差,使用雅可比矩阵代替一维导数,多变量方程的不动点迭代法,不动点存在唯一性定理(压缩映射原理),p阶收敛,牛顿迭代法等。

第8章 矩阵特征值计算

8.1 特征值性质和估计

特征值的基本性质,略。

格什戈林圆盘定理:设A = (a_{ij})_{n*n},则A的每一个特征值必属于下述某个圆盘之中:
| \lambda - a_{ii} | \leq r_i = \sum_{j=1,j\not= i}^n | a_{ij}|, \;\;\;\; i=1,2,\cdots,n
或者说,A的特征值都在复平面上n个圆盘的并集中。
如果Am个圆盘组成一个连通的并集S,且S与余下的n-m个圆盘是分离的,则S内恰包含Am个特征值。

圆盘定理可初步估计特征值的范围。

8.2 幂法及反幂法

幂法是一种计算矩阵主特征值(模最大的特征值)及其特征向量的迭代方法,适用于大型稀疏矩阵。
类似的,反幂法计算模最小的特征值及其特征向量,用于计算海森伯格矩阵或三对角阵的对于一个给定近似特征值的特征向量。

幂法:原理和编程都很简单,看书即可。前提:模最大的特征值只有一个(不是重特征值)。

加速方法:
1.原点平移法,引进矩阵B=A-pI,B是A的平移,特征值也对应平移且特征向量不变,可以改变模最大的特征值、收敛速度。缺点:p的确定困难。
2.瑞利商加速。

反幂法:计算A的模最小特征值,等价于计算A^{-1}的模最大特征值的倒数。(实际上不需要计算A^{-1},只有解线性方程组即可)
若已知某个特征值的近似值p(由其他方法求出),可以对矩阵平移B=A-pI,则B就有一个接近0的特征值,可使用反幂法求这个值和特征向量,也就是求得原矩阵A的更精确的特征值和对应的特征向量。

8.3 正交变换与矩阵分解

正交变换是计算矩阵特征值的有力工具,一般可以先通过正交变换将矩阵变为上三角矩阵或近似上三角的矩阵,再用迭代法快速求出所有特征值(QR方法)。
本节介绍豪斯霍尔德(Householder)变换吉文斯(Givens)变换

豪斯霍尔德变换:思路:使用反射,将x反射到e_1。计算量稍大,但一次改变向量的所有分量。

吉文斯变换:思路:使用旋转,将(\cdots,x_i,\cdots,x_j,\cdots)旋转到(\cdots,1,\cdots,0,\cdots)。计算量较小,但一次只改变两个分量。

矩阵的QR分解和舒尔分解
QR分解:设n阶方阵A非奇异,则可以通过一系列变换(比如n-1次豪斯霍尔德变换)将其变为上三角阵,即PA=R,其中P为正交矩阵,R为上三角矩阵。因此令Q=P^{-1},就有A=QR的分解。我们只需再求上三角矩阵R的特征值,就可以得到A的特征值。
性质:对n阶非奇异方阵A,QR分解存在且当R对角元为正时唯一。

实舒尔分解:解决实矩阵的复特征值问题,实矩阵可以约化为上海森伯格矩阵(上三角和对角线元素的下面一个元素可以非零,其余是零)。
Q^TAQ = R
R为类上三角矩阵,其对角块R_{ii}为一阶或二阶方阵。一阶对应实特征值,二阶对应两个共轭复特征值。计算Q的方法也是使用正交变换。

8.4 QR方法

QR方法是一种变换方法,是计算一般矩阵(中小型矩阵)全部特征值问题的最有效按方法之一。
对一般矩阵,先用豪斯霍尔德方法将A化为上海森伯格矩阵B,再用QR方法计算B的全部特征值(,再使用反幂法提高精度并求出对应的特征向量)。
迭代过程:对A_k做QR分解:A_k=Q_kR_k,再令A_{k+1}=R_kQ_K = Q_k^T R_kQ_k

具体细节,比如迭代法何时终止的判定,略。
还有在其基础上的改进,比如带原点位移的QR方法隐式QR方法,比较复杂,略。

第9章 常微分方程初值问题数值方法

9.1 引言

常微分方程初值问题:
y' = f(x,y),\;\;\;\; x\in [x_0,b] \\ y(x_0) = y_0
称方程的解y(x)为积分曲线。

李普希兹条件:如果存在实数L>0,使得
| f(x,y_1) - f(x,y_2) | \le L | y_1 -y_2|, \;\;\;\; \forall y_1,y_2 \in R
则称f关于y满足李普希兹条件。
这个条件用于常微分方程解的存在唯一性,具体参考《常微分方程》。在本章中也用于数值方法的收敛性证明。

对于常微分方程初值问题的数值方法,就是寻求一系列离散节点x_1<x_2<\cdots<x_n<\cdots上的近似值y_1<y_2<\cdots<y_n<\cdots
为了简便,我们取x_i = x_0 + i*h的等距节点,h称为步长。

一般使用前面的点(已知的或已求出近似值的)来计算后面的点。使用前一个点的叫做单步法,使用前n个点的叫做多步法,这和7.5节中的多点迭代法类似。
本章我们需要研究的问题有:计算方法(公式),公式的局部阶段误差和阶,误差估计及收敛性,稳定性。

9.2 简单的数值方法

理论上有y_{n+1} = y_n + \int_{x_n}^{x_{n+1}} y'(x)dx = y_n + \int_{x_n}^{x_{n+1}} f(x,y(x))dx

欧拉法y_{n+1} = y_n + hf(x_n,y_n),即使用左矩形公式计算数值积分,\int_{x_n}^{x_{n+1}} f(x,y(x))dx \approx hf(x_n,y_n)
这是最简单的思路最清晰的方法,计算量小,当然误差也比较大(截断误差相当于数值积分的误差)。
这是一个不需要解方程就可以直接计算的方法,即y_{n+1} = y_n + hf(x_n,y_n)等号右边没有未知量,这样的方法称为显式方法。

后退欧拉法(隐式欧拉法)y_{n+1} = y_n + hf(x_{n+1},y_{n+1}),思路与欧拉法相同,使用右矩形公式计算数值积分。此时等号右边有未知量,需要通过解方程来求得y_{n+1},这样的方法称为隐式方法,类似于隐函数和显函数的关系。
本方法的收敛性证明用到李普希兹条件
一般来说,隐式方法虽然计算量稍大,但稳定性好。

梯形方法:类似的,可以使用梯形公式计算数值积分。

改进欧拉公式:每次递推分成两部:先用欧拉公式求得“预测值”,再用梯形公式求得“校正值”,将梯形公式的隐式方法改为了两部显式计算,减少了计算量。

局部截断误差:加上目前的数据是正确的,下一个节点的数值解和精确值的差称为局部截断误差,以此来度量算法的收敛性和收敛速度。一般局部截断误差与步长h相关,即步长越小误差越小,若与h^{p+1}线性相关则称为p阶精度。

9.3 龙格-库塔方法

为了提高精度而增加了使用的点,增多的是[x_n,x_{n+1}]内的点,不是x_n之前的节点,因此还是单步法。类似改进的欧拉法,公式如下;
y_{n+1} = y_n + h\varphi(x_n,y_n,h) \\ \varphi(x_n,y_n,h) = \sum_{i=1}^r c_i K_i \\ K_1 = f(x_n,y_n) \\ K_i = f \left (x_n+\lambda_i h, y_n+h \sum_{j=1}^{i-1}\mu_{ij}K_j \right)
称为r阶显式龙格库塔法,简称RK法,c_i,\lambda_i,\mu_{ij}都是常参数。其中r=1时即为欧拉法,r=2时的一种为改进的欧拉法。

二阶显式R-K方法:以二阶为例说明参数的选择:希望确定这些参数使得p阶数尽量高,可使用待定系数法,使用泰勒公式并令几个低阶的误差项参数为零,从而得到几个方程,解得参数。由于方程的解不一定唯一,参数不唯一。
三阶、四阶类似。四阶显式龙格库塔方法比较常用。
由于龙格库塔方法的推导基于泰勒展开,因而它要求解具有较好的光滑性。

变步长的龙格库塔方法:略。

9.4 单步法的收敛性与稳定性

收敛性:步长趋于0时误差也趋于0,则称为具有收敛性。

定理:若单步法具有p阶精度且增量函数\varphi(x,y,h)关于y满足李普希兹条件,初值也准确,则y(x_n) - y = O(h^p),整体截断误差p阶收敛。
因此,一般通过判断是否符合李普希兹条件来判断是否收敛。

稳定性:若计算出现的微小误差不会恶性增长以至于“淹没”真解,则称具有稳定性。
而稳定性不仅与方法有关,有时也与h的取值范围有关(在《偏微分方程数值解》中体现比较明显,有的方法是全局收敛,有的方法是条件收敛,即必须步长/步长比较小时才收敛),收敛时h的取值范围叫做绝对稳定域/绝对稳定区间

9.5 线性多步法

思想:充分利用前面多步的节点值信息来预测下一个节点的值,期望会获得较高的精度。
方法:基于数值积分的方法或基于泰勒展开的方法。

线性多步法的一般公式:
y_{n+k} = \sum_{i=0}^{k-1} \alpha_i y_{n+i} +h\sum_{i=0}^k \beta_i f_{n+i}
类似龙格库塔公式,可以使用待定系数法和泰勒展开确定参数。

阿当姆斯显式与隐式公式:公式如下
y_{n+k} = y_{n+k-1} + h\sum_{i=0}^k\beta_if_{n+i}
其中\beta_k=0时为显式方法,反之为隐式方法。也称为阿当姆斯-巴什福斯(AB)公式阿当姆斯-蒙尔顿(AM)公式。可通过泰勒展开并求解线性方程组确定其系数。

还有其他的线性多步法公式,比如米尔尼方法辛普森方法汉明方法等。

预测校正方法:类似改进的欧拉方法,分为预测(显式公式)和校正(隐式公式)两步的多步方法,比如阿当姆斯四阶预测校正格式(PECE)

一般先用四阶龙格库塔方法(单步法里最常用的,效果较好但计算量大)算出开始的几个点,之后使用多步预测校正方法(计算量小)。

9.6 线性多步法的收敛性与稳定性

类似的,分为收敛性和稳定性,略。

9.7 一阶方程组与刚性方程组

略。

最后编辑于
©著作权归作者所有,转载或内容合作请联系作者
  • 序言:七十年代末,一起剥皮案震惊了整个滨河市,随后出现的几起案子,更是在滨河造成了极大的恐慌,老刑警刘岩,带你破解...
    沈念sama阅读 194,457评论 5 459
  • 序言:滨河连续发生了三起死亡事件,死亡现场离奇诡异,居然都是意外死亡,警方通过查阅死者的电脑和手机,发现死者居然都...
    沈念sama阅读 81,837评论 2 371
  • 文/潘晓璐 我一进店门,熙熙楼的掌柜王于贵愁眉苦脸地迎上来,“玉大人,你说我怎么就摊上这事。” “怎么了?”我有些...
    开封第一讲书人阅读 141,696评论 0 319
  • 文/不坏的土叔 我叫张陵,是天一观的道长。 经常有香客问我,道长,这世上最难降的妖魔是什么? 我笑而不...
    开封第一讲书人阅读 52,183评论 1 263
  • 正文 为了忘掉前任,我火速办了婚礼,结果婚礼上,老公的妹妹穿的比我还像新娘。我一直安慰自己,他们只是感情好,可当我...
    茶点故事阅读 61,057评论 4 355
  • 文/花漫 我一把揭开白布。 她就那样静静地躺着,像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上,一...
    开封第一讲书人阅读 46,105评论 1 272
  • 那天,我揣着相机与录音,去河边找鬼。 笑死,一个胖子当着我的面吹牛,可吹牛的内容都是我干的。 我是一名探鬼主播,决...
    沈念sama阅读 36,520评论 3 381
  • 文/苍兰香墨 我猛地睁开眼,长吁一口气:“原来是场噩梦啊……” “哼!你这毒妇竟也来了?” 一声冷哼从身侧响起,我...
    开封第一讲书人阅读 35,211评论 0 253
  • 序言:老挝万荣一对情侣失踪,失踪者是张志新(化名)和其女友刘颖,没想到半个月后,有当地人在树林里发现了一具尸体,经...
    沈念sama阅读 39,482评论 1 290
  • 正文 独居荒郊野岭守林人离奇死亡,尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 34,574评论 2 309
  • 正文 我和宋清朗相恋三年,在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 36,353评论 1 326
  • 序言:一个原本活蹦乱跳的男人离奇死亡,死状恐怖,灵堂内的尸体忽然破棺而出,到底是诈尸还是另有隐情,我是刑警宁泽,带...
    沈念sama阅读 32,213评论 3 312
  • 正文 年R本政府宣布,位于F岛的核电站,受9级特大地震影响,放射性物质发生泄漏。R本人自食恶果不足惜,却给世界环境...
    茶点故事阅读 37,576评论 3 298
  • 文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹,春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 28,897评论 0 17
  • 文/苍兰香墨 我抬头看了看天上的太阳。三九已至,却和暖如春,着一层夹袄步出监牢的瞬间,已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 30,174评论 1 250
  • 我被黑心中介骗来泰国打工, 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留,地道东北人。 一个月前我还...
    沈念sama阅读 41,489评论 2 341
  • 正文 我出身青楼,却偏偏与公主长得像,于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子,可洞房花烛夜当晚...
    茶点故事阅读 40,683评论 2 335