地球物理计算方法第一章

第一章

线性方程组的数值解法2/47本章内容§1.1

引言§1.2

高斯消去法§1.3

选主元素的高斯消去法§1.4

矩阵的三角分解§1.5解三对角线方程组的追赶法§1.6解对称正定矩阵方程组的平方根法§1.7向量和矩阵的范数§1.8解线性方程组的迭代法§1.9病态方程组和迭代改善法3/47本章要求1.理解高斯消去法的基本原理，熟练掌握高斯主元消去法；2.理解矩阵的三角分解；3.掌握解三对角线方程组的追赶法，掌握平方根法；4.了解矩阵范数、条件数。5.熟悉简单迭代法及其收敛条件的使用;6.熟悉Jacobi迭代法及其相应的Seidel迭代法的计算公式以及它们的收敛条件;7.熟悉SOR方法的计算公式及其收敛条件。4/47§1.1引言本节内容一.线性方程组解法二.直接法与迭代法比较5/47§1.1引言一.线性方程组解法工程中几乎有一半的问题涉及到线性方程组的求解设n阶线性方程组6/47§1.1引言A

称为方程组的系数矩阵，当是n阶非奇异矩阵时，既A0，此时方程组有唯一解。

X阵是解向量，B阵是常向量。在线性代数中学过用克莱姆法则求解，它是一种直接法(属于解析法)，但随着n计算工作量7/47§1.1引言自学并理解8/47§1.1引言二.直接法与迭代法比较各有优缺点解方程组的数值解法：直接法与迭代法优缺点比较直接法：经有限次计算得准确解(在无舍入误差下),实际上舍入误差客观存在,得到的依然还是近似解。由于受到计算机存储容量的限制，一般来说，仅适于系数矩阵阶数不太高的问题，其工作量(计算量)较小、精度较高，但程序设计复杂。迭代法将问题构成一个无穷序列，逼近准确解。主要用于某些系数矩阵阶数较高的问题，一般来说，程序较为简单、易于编程，但存在收敛性及收敛速度的问题，只对具有某些性质的系数矩阵的方程组才适用。工作量有时较大。9/47§1.1引言实际计算时，应根据问题的特点和要求来决定方法的取舍。本章介绍的求解线性代数方程组的直接法有Gauss(高斯)消元法和LU分解等；迭代法有Jacobi(雅可比)迭代和Gauss-Seidel(高斯-赛德尔)迭代。10/47§1.2高斯消去法本节内容一.引言二.例子三.顺序Gauss消去法四.Gauss消去法计算量返回章节目录顺序高斯消去法11/47一.引言是一个古老的求解线性方程组的方法。改进和变形得到高斯选主元素消去法（第3节）、三角分解法（第4节）n元线性方程组的直接解法。§1.2高斯消去法（式1）12/47方程组(1)的矩阵形式为Ax=b

其中§1.2高斯消去法13/47线性代数：方法不好时工作量非常大，工作量小的方法是Gauss消去法。

Cramer法则是一种不实用的直接法，本章介绍几种实用的直接法。

Gauss消去法是一种规则化的加减消元法，其基本思想是通过逐次消元计算，把一般线性方程组的求解转化为等价的上三角形方程组的求解。§1.2高斯消去法14/47二.例子为清楚起见，先看一简单例子。考虑线性方程组1.消去后两个方程中的x1得：2.再消去最后一个方程的x2得：3.消元结束,经过回代得解:§1.2高斯消去法15/47上述求解的消元过程可用矩阵表示为:

这是高斯消去法的计算形式,新的增广矩阵对应的线性方程组就是上三角形方程组,可进行回代求解§1.2高斯消去法16/47三.

求解线性方程组(1)的顺序Gauss消去法记则,线性方程组(1)的增广矩阵为§1.2高斯消去法17/47§1.2高斯消去法

18/47§1.2高斯消去法

19/47§1.2高斯消去法

20/47（式2）§1.2高斯消去法

（式3）21/47§1.2高斯消去法22/47结论：

（1）高斯消去法分消元、回代两过程。

（2）从矩阵分解角度看：

消去是解一个下三角方程组，

回代是解一个上三角方程组。

（3）消去法顺利进行必须满足

akk

(k)0,(k=1,2,…,n)，若出现akk

(k)=0，则可交换行列后再进行消元。§1.2高斯消去法23/47§1.2高斯消去法四.Gauss消去法计算量(n-k)(n-k)224/47§1.2高斯消去法25/47§1.2高斯消去法乘除法耗时大大多于加减法耗时，故高斯消元法的计算量为O(n3)。n=20时,顺序Gauss消去法只需3060次乘除法运算。顺序Gauss消去法通常也简称为Gauss消去法。26/47§1.2高斯消去法27/47§1.3

选主元素的高斯消去法本节内容一.问题提出二.选主元素消去法三.高斯列主元消去法N-S图四.高斯-约当消去法返回章节目录28/47一.问题提出§1.3

选主元素的高斯消去法29/47§1.3

选主元素的高斯消去法例1：单精度解方程组精确解为和8个8个用Gauss消去法计算：8个小主元

/*Smallpivotelement*/

可能导致计算失败。30/47§1.3

选主元素的高斯消去法

例2：解线性方程组（用4位十进制浮点计算）用Cramer法则可得精度较高的解（精确解）

x1*=1.00010，x2*=0.9999031/47

解用顺序Gauss消去法,消元得回代得解：x1=0.00，x2=1.00与精确解相比，该结果相当糟糕原因是：在求行乘数时用了很小的数0.0001作除数主元9999§1.3

选主元素的高斯消去法32/47用顺序Gauss消去法,消元得回代得解：x1=1.00，x2=1.00若将方程组改写成（将1,2行交换）:这是一个相当不错的结果0.9999§1.3

选主元素的高斯消去法33/47例3：解线性方程组(用8位十进制尾数的浮点数计算)解：这个方程组和例2一样，若用Gauss消去法计算会有小数作除数的现象，若采用换行的技巧，则可避免§1.3

选主元素的高斯消去法34/47绝对值最大,不需换行§1.3

选主元素的高斯消去法35/47§1.3

选主元素的高斯消去法36/47§1.3

选主元素的高斯消去法二.选主元素消去法37/47§1.3

选主元素的高斯消去法38/47

给定线性方程组Ax=b，记A(1)=A，b(1)=b，列主元Gauss消去法的具体过程如下：首先在增广矩阵B(1)=(A(1),b(1))的第一列元素中，取然后进行第一步消元得增广矩阵B(2)=(A(2),b(2))。再在矩阵B(2)=(A(2),b(2))的第二列元素中，取然后进行第二步消元得增广矩阵B(3)=(A(3),b(3))。按此方法继续进行下去,经过n-1步选主元和消元运算,得到增广矩阵B(n)=(A(n),b(n)).则方程组A(n)x=b(n)是与原方程组等价的上三角形方程组,可进行回代求解.

易证，只要|A|0,列主元Gauss消去法就可顺利进行。§1.3

选主元素的高斯消去法39/47

方程组具有四位有效数字的精确解为

x1*=17.46，x2*=-45.76，x3*=5.546例4采用4位十进制浮点计算,分别用顺序Gauss消去法和列主元Gauss消去法求解线性方程组:§1.3

选主元素的高斯消去法40/47

解1.用顺序Gauss消去法求解，消元过程为回代得：x3=5.546，x2=100.0，x1=-104.0§1.3

选主元素的高斯消去法41/472.用列主元Gauss消去法求解，消元过程为§1.3

选主元素的高斯消去法42/47回代得：x3=5.545，x2=-45.77，x1=17.46§1.3

选主元素的高斯消去法43/47

列主元Gauss消去法是在每一步消元前，在主元所在的一列选取绝对值最大的元素作为主元素。而全主元Gauss消去法是在每一步消元前，在所有元素中选取绝对值最大的元素作为主元素。但由于运算量大增，实际应用中并不经常使用。§1.3

选主元素的高斯消去法44/47三.高斯列主元消去法N-S图消元次数K=1k<=n?选列主元amk|amk|<epsNrmrk进行第k步消元输出方程组的解X建立增广矩阵A(n,n+1)m=kYNYk=k+1回代求解停止§1.3

选主元素的高斯消去法45/47四.高斯-约当消去法（Gauss-Jordan）这是一种无回代的方法，使计算工作量减小，在第k次消元时不仅把akk(k)

化成1，不仅将aik（i=k+1，……，n）化成0，而且也将aik（i=1，2，……，k-1）也化为0，最后消元过程结束时的系数矩阵为得消元公式

akj

(k)=akj(k-1)/akk(k-1)

（j=1,2,…,n+1）

aij

(k)=aij

(k)-aik(k-1)akj(k)

（ik,i=1,2,…,n，j=1,2,…,n+1）自学了解§1.3

选主元素的高斯消去法46/47§1.4

矩阵的三角分解本节内容一.Gauss消去法的矩阵表示二.Doolittle分解法返回章节目录47/47§1.4

矩阵的三角分解一.Gauss消去法的矩阵表示对矩阵48/47§1.4

矩阵的三角分解49/47§1.4

矩阵的三角分解50/47§1.4

矩阵的三角分解51/47也就是：

A(n)=Ln-1A(n-1)=Ln-1Ln-2A(n-2)=…=Ln-1Ln-2…L2L1A(1)

其中：所以有：A=A(1)=L1-1L2-1…Ln-1-1A(n)=LU

其中：L=L1-1L2-1…Ln-1-1,U=A(n)§1.4

矩阵的三角分解52/47L称为单位下三角矩阵;U是上三角矩阵.式A=LU称为矩阵A的三角分解。§1.4

矩阵的三角分解53/47二、Doolittle分解法定理2：设n阶方阵A的各阶顺序主子式不为零,则存在唯一单位下三角矩阵L和上三角矩阵U使A=LU。证明只证唯一性,设有两种分解则有所以得于是Ax=bLUx=b

令Ux=y得§1.4

矩阵的三角分解54/47K=1§1.4

矩阵的三角分解55/47对K=2,3,…,n计算（1）求U的第i行元素§1.4

矩阵的三角分解56/47对K=2,3,…,n计算（2）求L的第i列元素§1.4

矩阵的三角分解57/47§1.4

矩阵的三角分解58/47§1.4

矩阵的三角分解59/47这就是求解方程组Ax=b的Doolittle三角分解方法。§1.4

矩阵的三角分解60/47例利用三角分解方法解线性方程组

解因为§1.4

矩阵的三角分解61/47§1.4

矩阵的三角分解先解,得再解,得62/47

解线性方程组Ax=b

的Doolittle

三角分解法的计算量约为1/3n3,与Gauss消去法基本相同。其优点在于求一系列同系数的线性方程组Ax=bk,(k=1,2,…,m)时,可大大节省运算量。例如,求上例中矩阵A的逆矩阵。可分别取常向量

b1=(1,0,0)T,b2=(0,1,0)T,b3=(0,0,1)T

§1.4

矩阵的三角分解63/47§1.4

矩阵的三角分解64/47

为了提高数值稳定性,可考虑列主元三角分解法,设已完成A=LU的k-1步分解计算,矩阵分解成§1.4

矩阵的三角分解65/47§1.4

矩阵的三角分解66/47§6.4

矩阵的三角分解例如，用列主元三角分解解上例中方程组。则有67/47§1.5

解三对角线方程组的追赶法本节内容一.三对角线矩阵二.Crount分解三.例返回章节目录68/47一.

有一类方程组，在插值问题和边值问题中有着重要的作用，即三对角线方程组，其形式为：§1.5

解三对角线方程组的追赶法对角元前次对角线后次对角线69/47§1.5

解三对角线方程组的追赶法70/47§1.5

解三对角线方程组的追赶法前面已经讲过71/47二对角矩阵§1.5

解三对角线方程组的追赶法注意：和教材上字母不同72/47§1.5

解三对角线方程组的追赶法73/47§1.5

解三对角线方程组的追赶法74/47追，相当于消元过程§1.5

解三对角线方程组的追赶法75/47式1—式3叫追赶法，也称Thomas法。工作量小，非常有效。赶，相当于回代过程§1.5

解三对角线方程组的追赶法76/47§5.7

解三对角线方程组的追赶法77/47§5.7

解三对角线方程组的追赶法78/47当满足条件

|a1|>|c1|>0；

|an|>|dn|>0；

|ai||ci|+|di|，cidi0，i=2,3,…,n-1。时，追赶法是数值稳定的追赶法具有：计算程序简单，存贮量少，计算量小的优点。§1.5

解三对角线方程组的追赶法79/47§1.6解对称正定矩阵方程组的平方根法本节内容一.对称正定矩阵二.平方根法——乔来斯金分解三.LDLT分解返回章节目录80/47§1.6解对称正定矩阵方程组的平方根法81/47§1.6解对称正定矩阵方程组的平方根法

对称正定阵的几个重要性质(2)A的顺序主子阵Ak

亦对称正定(3)A的特征值i>0

(1)A1

亦对称正定，且aii>0(4)A的全部顺序主子式

det(Ak)>082/47§1.6解对称正定矩阵方程组的平方根法

83/47§1.6解对称正定矩阵方程组的平方根法84/47§1.6解对称正定矩阵方程组的平方根法85/47§1.6解对称正定矩阵方程组的平方根法86/47§1.6解对称正定矩阵方程组的平方根法87/47§1.6解对称正定矩阵方程组的平方根法88/47§6.6解对称正定矩阵方程组的平方根法√44/2=289/47§6.6解对称正定矩阵方程组的平方根法

平方根法是求对称正定系数线性方程组的三角分解法，对称正定矩阵的Cholesky

分解的计算量和存贮量均约为一般矩阵的LU

分解的一半。且Cholesky

分解具有数值稳定性。90/47§1.6解对称正定矩阵方程组的平方根法

91/47§1.6解对称正定矩阵方程组的平方根法92/47§1.7向量和矩阵的范数本节内容一.向量范数二.矩阵范数三.矩阵A的谱半径返回章节目录93/47§1.7向量和矩阵的范数

94/47§1.7向量和矩阵的范数95/47§1.7向量和矩阵的范数96/47§1.7向量和矩阵的范数97/47§1.7向量和矩阵的范数98/47§1.7向量和矩阵的范数

99/47§1.7向量和矩阵的范数100/47§1.7向量和矩阵的范数101/47§1.7向量和矩阵的范数102/47§6.7向量和矩阵的范数103/47§1.7向量和矩阵的范数104/47§1.7向量和矩阵的范数105/47§1.8解线性方程组的迭代法本节内容一.简单迭代法二.雅可比(Jacobi)迭代法三.塞德尔(Seidel)迭代法四.

逐次超松弛(SOR)迭代法五.

例子返回章节目录解大型稀疏线性方程组106/47§1.8解线性方程组的迭代法一.简单迭代法

107/47§1.8解线性方程组的迭代法一阶定常迭代法108/47

特征值谱半径§1.8解线性方程组的迭代法

109/47§6.8解线性方程组的迭代法有兴趣同学自学110/47二.雅可比(Jacobi)迭代法1.迭代法§1.8解线性方程组的迭代法111/47§1.8解线性方程组的迭代法112/47§1.8解线性方程组的迭代法113/47§1.8解线性方程组的迭代法114/47§1.8解线性方程组的迭代法115/47§1.8解线性方程组的迭代法2.收敛性范数116/47§1.8解线性方程组的迭代法117/47§1.8解线性方程组的迭代法三.塞德尔(Seidel)迭代法1.迭代法1118/47§1.8解线性方程组的迭代法119/47§1.8解线性方程组的迭代法2.收敛性1120/47§1.8解线性方程组的迭代法3.迭代法2121/47§1.8解线性方程组的迭代法122/47§1.8解线性方程组的迭代法4.收敛性2123/47四.

逐次超松弛(SOR)迭代法1.迭代法§1.8解线性方程组的迭代法124/47§1.8解线性方程组的迭代法2.收敛性125/47§1.8解线性方程组的迭代法五.

例子例1126/47§1.8解线性方程组的迭代法127/47§6.8解线性方程组的迭代法可见，迭代序列逐次收敛于方程组的解，而且迭代7次得到精确到小数点后两位的近似解。kx1(k)x2(k)x3(k)‖x(k)-x*‖0123456701.41.110.9290.99061.011591.0002510.998236400.51.201.0550.96450.99531.0057951.000125501.41.110.9290.99061.011591.0002510.998236410.50.20.0710.03550.011590.0057950.0017636计算结果列表如下:128/47G-S迭代法的计算公式为:§1.8解线性方程组的迭代法129/47同样取初始向量x(0)=(0,0,0)T，计算结果为：

由计算结果可见，G-S迭代法收敛较快.取精确到小数点后两位的近似解，G-S迭代法只需迭代3次，而J迭代法需要迭代7次。kx1(k)x2(k)x3(k)‖x(k)-x*‖012301.41.06340.995104400.781.020480.9952756801.0260.9875161.0019068610.40.06340.0048956§1.8解线性方程组的迭代法130/47若使‖x(k)–x*‖<,只需,即可以事先估计达到某一精度需要迭代多少步。§1.8解线性方程组的迭代法131/47例2用J迭代法求例1中方程组的解，取x(0)=(0,0,0)T，若使误差x(k)-x*<10-5，问需要迭代多少次？解由例1知，x(1)=(1.4,0.5,1.4)T，于是有，x(1)-x(0)=1.4，B=0.5。k应满足故取k=19，即需要迭代1