数值计算方法课件-第五章-线性方程组的数值解法

1Chapter5NumericalSolutionofLinearEquations线性方程组的数值解法23直接法:经过有限步运算后可求得方程组精确解的方法(不计舍入误差!)迭代法：从解的某个近似值出发，通过构造一个无穷序列去逼近精确解的方法。分为两类：

逐次逼近法(一般有限步内得不到精确解)

共轭斜量法（不考虑计算过程的舍入误差，只用有限步就收敛于方程组的精确解）解线性方程组的两类方法

4

DirectMethodforSolvingLinearSystems直接解法：它是一类精确方法，即若不考虑计算过程中的舍入误差，那么通过有限步运算可以获得方程解的精确结果.

Gauss逐步（顺序）消去法、Gauss主元素法、矩阵分解法等；5

IterativeMethodsforSolvingLinearSystems

迭代解法：所谓迭代方法，就是构造某种极限过程去逐步逼近方程组的解.

经典迭代法有：

迭代法、迭代法、逐次超松弛（SOR）迭代法等；6§5.1解线性方程组的直接法

(DirectMethodforSolvingLinearSystems)高斯消去法

(GaussianElimination)思路首先将A化为上三角阵(upper-triangularmatrix)，此过程称为消去过程，再求解如下形状的方程组，此过程称为回代求解

(backwardsubstitution)。=5.1.1求解的高斯消去法和选主元高斯消去法7将增广矩阵的第i行+li1

第1行，得到：消去过程：第一步：设，计算因子其中8将增广矩阵的第i行+lik

第k行，得到：其中第k步：设，计算因子9定理：若A的所有顺序主子式均不为0，则高斯消去法能顺序进行消元，得到唯一解。回代过程：共进行n

1步，得到10运算量

（AmountofComputation）（1）用克莱姆（Cramer）法则求解n阶线性方程组

每个行列式由n!项相加，而每项包含了n个因子相乘，乘法运算次数为(n-1)n!次.仅考虑乘(除)法运算,计算解向量包括计算n+1个行列式和n次除法运算,乘(除)法运算次数N=(n+1)(n-1)n!+n.当n=8时,N＝200，000011（2）高斯消去法:

第1个消去步，计算li1(i=2,3，…，n)，有n-1次除法运算.使aij(1)变为aij(2)

以及使bi(1)变为bi(2)有n(n-1)次乘法运算.

第k个消去步，有n-k次除法运算、(n-k+1)(n-k)次乘法运算.乘法运算总次数为：除法运算总次数为:(n-1)+…+1=n(n-1)/212回代过程的计算除法运算次数为n次.乘法运算的总次数为

n+(n-1)+…+1=n(n-1)/2次Gauss消去法除法运算次数为:n(n-1)/2+n=n(n+1)/2，乘法运算次数为:n(n-1)(n+1)/3+n(n-1)/2=n(n-1)(2n+5)/6，通常也说Gauss消去法的运算次数与n3同阶，记为O(n3)13二、选主元消去法在高斯消去法消去过程中可能出现的情况，这时高斯消去法将无法进行；即使主元素但很小，其作除数，也会导致其它元素数量级的严重增长和舍入误差的扩散.14Ex：单精度解方程组

精确解为和8个8个用Gauss消去法计算：8个小主元

Smallpivotelement15列主元消去法在第k步消元前，在系数矩阵第k列的对角线以下的元素中找出绝对值最大的元。

若p≠k，交换第k个与第p个方程后，再继续消去计算.

这种方法称为列主元Gauss消去法。

列主元Gauss消去法保证了｜lik｜≤1(i=k+1,k+2,…，n).为避免这种情况的发生,可通过交换方程的次序，选取绝对值大的元素作主元.基于这种思想导出了“选主元消去法”16

全主元消去法在第k步消去前，在系数矩阵右下角的n-k+1阶主子阵中，选绝对值最大的元素作为主元素。(1)Ifpk

then交换第k行与第p行;Ifqk

then交换第k列与第q列;(2)

消元注：列交换改变了xi

的顺序，须记录交换次序，解完后再换回来。175.1.2

三角分解法

高斯消元法的矩阵形式

每一步消去过程相当于左乘初等变换矩阵Lk1819A

的

LU

分解(LUfactorization)定理1.1：如果Gauss消去法能顺序进行消去，则矩阵A可进行三角分解，即A＝LU20注：（1）L

为单位下三角阵而U

为一般上三角阵的分解称为Doolittle

分解（2）L

为一般下三角阵而U

为单位上三角阵的分解称为Crout分解。

21Doolittle分解法：通过比较法直接导出L和

U的计算公式。思路22一般计算公式23LU分解求解线性方程组242526定理1.2：Gauss消去法能顺序进行消去的充要条件是矩阵A的所有顺序主子阵Ai

非奇异.27定理1.4：若A的所有顺序主子式均不为0，则A

的

LU

分解唯一（其中L

为单位下三角阵）。证明：由定理1.1和定理1.2可知，LU分解存在。下面证明唯一性。若不唯一，则可设A=L1U1=L2U2

，推出上三角矩阵对角线上为1的下三角矩阵285.1.4求解正定方程组的Cholesky方法对称正定阵A的几个重要性质：（1）A1

亦对称正定，且aii>0（2）A

的顺序主子阵

Ak

亦对称正定（3）A

的特征值i

>0

（4）A

的全部顺序主子式

det(Ak

)>029定理1.5设矩阵A对称正定，则存在唯一的对角元全为正的下三角阵G

使得

A＝GGT计算格式为305.1.5求解三对角方程组的追赶法定理：若A

为对角占优的三对角阵，且满足则方程组有唯一的LU分解。31直接比较等式两边的元素，可得到计算公式.第二步:追—即解：第三步:赶—即解：第一步:对A作Crout分解32

Def2:设｛xk｝是Rn上的向量序列，令xk=(xk1,xk2,…，xkn)Ｔ,k=1,2，….,

又设x*＝（x1*，x2*，…，xn*)Ｔ是Rn上的向量.

如果limxki=xi对所有的i=1,2,…，n成立，

那么,称向量x*是向量序列｛xk｝的极限，若一个向量序列有极限,称这个向量序列是收敛的.

对任意一种向量范数‖·‖而言，向量序列｛xk｝收敛于向量x*的充分必要条件是Theorem33

矩阵范数(matrixnorms)Def3:对任意,称||·||为Rmn空间的矩阵范数,

指||·||满足(1)-(3)：对任意(4)

||AB||||A||·||B||若还满足(4),称为相容的矩阵范数34Ex:设A＝(aij)∈M.定义证明：这样定义的非负实数不是相容的矩阵范数.证明：设从而35相容性（1）矩阵范数与矩阵范数的相容:‖AB‖≤‖A‖‖B‖（2）矩阵范数与向量范数相容设A∈M，‖A‖是矩阵范数,x∈Rn，‖x‖是向量范数.如果满足不等式:‖Ax‖≤‖A‖‖x‖则称矩阵范数‖A‖与向量范数‖x‖相容.36常用的算子范数：

由向量范数

||·||p

导出关于矩阵

A

Rnn

的p范数:则（行和范数）（列和范数）（谱范数(spectralnorm)

）对方阵和有:

(向量||·||2的直接推广)Frobenius范数:(operatornorm),又称为从属的矩阵范数:算子范数37Theorem对任意算子范数||·||有:证明：由算子范数的相容性，得到将任意一个特征根所对应的特征向量代入Theorem若A对称，则有:证明：若是A

的一个特征根，则2

必是A2

的特征根。又：对称矩阵的特征根为实数，即2(A)为非负实数，故得证。对某个

A

的特征根成立38Theorem若矩阵A

对某个算子范数满足||A||<1，则必有①.可逆；②.证明：①若不然，则有非零解，即存在非零向量使得②39

ErrorAnalysisforLinearsystemofEquations线性方程组的性态（误差分析）

思考:求解时,A

和的误差对解

有何影响？设A

精确，有误差，得到的解为，即绝对误差放大因子又相对误差放大因子40设精确，A有误差，得到的解为，即(只要A充分小，使得

是关键的误差放大因子，称为A的条件数，记为cond(A)

41注:cond(A)的具体大小与||·||的取法有关。

cond(A)取决于A，与解题方法无关。常用条件数有：cond(A)1cond(A)cond(A)2特别地，若A对称，则条件数的性质：

A可逆，则cond(A)p

1；

A可逆，R

则cond(

A)

=cond(A)

；

A正交，则cond(A)2=1；

A可逆，R正交，则cond(RA)2

=cond(AR)2

=cond(A)2

。42精确解为例计算cond

A)2。A1=解：考察A

的特征根39206>>1

测试病态程度：给一个扰动，其相对误差为此时精确解为2.0102>200%43例：Hilbert阵cond(H2)=27cond(H3)748cond(H6)=2.9106cond(Hn)asn注：一般判断矩阵是否病态，并不计算A1，而由经验得出。

行列式很大或很小（如某些行、列近似相关）；

元素间相差大数量级，且无规则；

主元消去过程中出现小主元；

特征值相差大数量级。44

对于大型线性方程组，Gauss消去法的计算工作量很大，用向前误差分析方法非常困难。下面介绍利用向后误差分析方法：Gauss消去法的舍入误差将实际计算过程的误差转换为原始数据的误差.45具体提法是：对两种扰动：考查对的解的影响，即的大小可以估计在十进制且16位字长的计算机上：46其中表示Gauss消元过程中的元素其中表示Gauss消元过程中的元素47特别，对Gauss消去法，只要不是很大，则Gauss消去法是数值稳定的结论：条件数越大，扰动对解的影响越大.48迭代法的基本思想Jacobi迭代和Gauss－Seidel迭代迭代法的收敛性超松弛迭代SOR§5.5解线性方程组的迭代法

49

迭代法的基本思想例：求解方程组其精确解是x*=(3,2,1)T。现将原方程组改写为简写为x=B0x+f，其中50任取初始值，如取x(0)=(0,0,0)T，代入x=B0x+f右边，若等式成立则求得方程组的解，否则，得新值x(1)=(x1(1),x2(1),x3(1))T=(2.5,3,3)T，再将x(1)代入x=B0x+f右边，反复计算，得一向量序列{x(k)}和迭代公式：简写为x(k+1)=B0x(k)+fk=0,1,2,……x(10)=(3.000032,1.999838,0.999813)T迭代到第10次时有||ε(10)||

∞=||x(10)–x*||=0.00018751Defination：（1）对于给定方程组x=Bx+f，用迭代公式x(k+1)=Bx(k)+f(k=0,1,2,……)逐步代入求近似解的方法称迭代法；（2）若k∞时limx(k)存在（记为x*），称此迭代法收敛，显然x*就是方程组的解，否则称迭代法发散；（3）B称为迭代矩阵.问题：

如何建立迭代格式？

收敛速度？

向量序列的收敛条件？

误差估计？52

求解思路

将等价改写为形式，建立迭代.从初值出发，得到序列.计算精度可控，特别适用于求解系数为大型稀疏矩阵的方程组。研究内容：

如何建立迭代格式？

收敛速度？

向量序列的收敛条件？

误差估计？53的收敛条件充分条件:||B||<1必要条件:？定义设：AAkk=lim是指ijkijkaa=)(lim对所有1i,jn

成立。迭代法的收敛性54对任意非零向量成立定理设存在唯一解，则从任意出发，迭代收敛0kB证明：Bk

0||Bk||0对任意非零向量成立从任意出发，记，则ask

收敛55定理：迭代格式收敛的充要条件为迭代矩阵的谱半径(B)<1.证:对任何n阶矩阵B，都存在非奇矩阵P，使

B=P–1JP其中，J为B的Jordan标准型。其中,Ji

为Jordan块56其中,λi

是矩阵B的特征值,由B=P–1JPBk=(P–1JP)(P–1JP)···(P–1JP)=P–1JkP迭代法x(k+1)=Bx(k)+f

收敛<=>(i=1,2,···,r)(i=1,2,···,r)谱半径(B)<157推论

若存在一个矩阵范数使得||B||=q<1,

则迭代收敛，且有下列误差估计：①②证明：①②58设Ax=b，A非奇异，且对角元不为零，将原方程组改写为Jacobi迭代法59又代入，反复继续，得迭代格式：称Jacobi迭代法.选取初始向量代入上面方程组右端得60Jacobi迭代法的矩阵表示：61故计算公式为：(i=1,2,……,n)，(k=0,1,2,……表迭代次数)矩阵表示：62则BJ=I-D-1

A =D-1(L+U)，

fJ=D-1b，称BJ为Jacobi迭代矩阵.（aii≠0）将方程组Ax＝b的系数矩阵A分解为：A=D-L-U63例1用Jacobi迭代法求解方程组,误差不超过10-4.解:646566依此类推,得方程组满足精度的解为x12迭代次数：12次

x4=3.02411.94780.9205d=0.1573x5=3.00031.98401.0010d=0.0914x6=2.99382.00001.0038d=0.0175x7=2.99902.00261.0031d=0.0059x8=3.00022.00060.9998d=0.0040x9=3.00031.99990.9997d=7.3612e-004x10=3.00001.99990.9999d=2.8918e-004x11=3.00002.00001.0000d=1.7669e-004x12=3.00002.00001.0000d=3.0647e-00567

Jacobi法Jacobi迭代方法写成矩阵形式：A=LUDBJacobi迭代阵68

算法:Jacobi迭代方法

求解给定初始值.Input:

方程和未知量的个数n;矩阵元素a[][];

元素b[];初始近似值X0[];误差容限TOL;

最大迭代次数Nmax.Output:

近似解X[]或失败的信息.Step1Setk=1;Step2While(kNmax)dosteps3-6

Step3Fori=1,…,n

Set;/*计算xk*/

Step4IfthenOutput(X[]);STOP;/*成功*/

Step5Fori=1,…,nSetX0[]=X[];/*更新X0*/

Step6Setk++;Step7Output(超过最大迭代次数);STOP./*失败*/69若在迭代时尽量利用最新信息，则可将迭代格式变为Gauss－Seidel迭代法称Gauss－Seidel迭代法.70计算公式：（i=2,3,……,n-1）(i=1,2,…,n)即71其中BG-S=(D-L)-1

U称Gauss－Seidel迭代矩阵.(D–L)x(k+1)=b＋

Ux(k)故Gauss－Seidel迭代格式：72例2用Jacobi迭代法求解方程组,误差不超过10-4.解:7374x1=2.50002.09091.2273d=3.4825x2=2.97732.02891.0041d=0.5305x3=3.00981.99680.9959d=0.0465x4=2.99981.99971.0002d=0.0112x5=2.99982.00011.0001d=3.9735e-004x6=3.00002.00001.0000d=1.9555e-004x7=3.00002.00001.0000d=1.1576e-005通过迭代,至第7步得到满足精度的解x7：从例1和例2可以看出,Gauss-Seidel迭代法的收敛速度比Jacobi迭代法要快.Jacobi迭代法和Gauss-Seidel迭代法统称为简单迭代法.75

Gauss-Seidel迭代方法…………写成矩阵形式：BGauss-Seidel

迭代阵76

定理

（充分条件）若A

为严格对角占优阵，则解的Jacobi和Gauss-Seidel迭代均收敛。证明：首先需要一个引理若A

为严格对角占优阵，则det(A)0，且所有的aii0。证明：若不然，即det(A)=0，则A

是奇异阵。存在非零向量使得记我们需要对Jacobi迭代和Gauss-Seidel迭代分别证明：任何一个||1

都不可能是对应迭代阵的特征根，即|IB|0

。Jacobi:BJ=D1(L+U)aii0如果||1则|IB|077

注：二种方法都存在收敛性问题。有例子表明：Gauss-Seidel法收敛时，Jacobi法可能不收敛；而Jacobi法收敛时，Gauss-Seidel法也可能不收敛。78例：给出方程组其中问:Jacobi迭代法和Gauss-Seidel迭代法是否收敛？解：对79而即所以，对Jacobi方法收敛，G-S方法发散.

同理，对于其中80即得而81则所以，对Jacobi方法发散，G-S方法收敛.说明，Jacobi方法和G-S方法的收敛条件大多数是互不包含的.82例2判别下列方程组用Jacobi法和Gauss-Seidel法求解是否收敛:解:(1)求Jacobi法的迭代矩阵83所以即Jaobi迭代法收敛。(2)求Gauss-Seidel法的迭代矩阵：显然BJ的几种常用算子范数||BJ||>1，故用其特征值判断。84所以Gauss-Seidel迭代法发散。注1：Gauss-Seidel迭代法收敛速度比Jacobi迭代法要高，但例2却说明Gauss-Seidel迭代法发散时而Jacobi迭代法却收敛，因此,不能说Gauss-Seidel迭代法比Jacobi迭代法更好.可得故85注2（2）Jacobi迭代法和Gauss-Seidel迭代法的收敛性没有必然的联系：即当Gauss-Seidel法收敛时，Jacobi法可能不收敛；而Jacobi法收敛时，Gauss-Seidel法也可能不收敛.（1）Jacobi迭代法和Gauss-Seidel迭代法的迭代矩阵不同：BJ=D-1(L+U),BG-S=(D-L)-1U86定义设A=(aij)n×nRn×n

，若

(i=1,2,…,n)则称A为对角占优矩阵，若不等式严格成立，则称A为严格对角占优矩阵.迭代法收敛的其他结论：定理若Ax=b中A为严格对角占优矩阵，则Jacobi迭代和Gauss－Seidel迭代均收敛。证明：因为系数矩阵A严格对角占优，所以87故Jacobi迭代法收敛(1)对于Jacobi迭代法,其迭代矩阵为88即从而因此由于可得(2)对于G—S迭代法,其迭代矩阵为89矛盾由前述定理知，G—S迭代法收敛定理

若A为对称正定阵，则Gauss－Seidel迭代收敛。90定理

若线性方程组的系数矩阵是对角元素为正的实对称阵，则Jacobi方法收敛的充分必要条件是和同时正定.证明.如果A是有正对角元的对称矩阵，aii>0记D=diag(a11，a22,…，ann).

91对于BJ=I-D-1A,有(1)若Jacobi迭代法收敛，则

(BJ)<1，即

即D－1/2AD-1/2是正定矩阵，从而A也是正定矩阵.现考虑矩阵2D－A即2D－A与的特征值的符号相同92的特征值为2－i因此2D-A的特征值为正，从而2D-A为正定矩阵充分性：A为正定矩阵I－D－1A的特征值全小于12D-A为正定矩阵

I－D－1A的特征值全大于－1（I－D－1A）<1（BJ）<1

93例：给定，其中⑴证明当时,对称正定(G-S迭代方法收敛).⑵证明当时,Jacobi迭代方法收敛，否则发散.94（1）为对称阵证明：时,对称正定，从而G-S迭代方法收敛.95（2）由定理1，是对角元素为正的实对称阵，和同时正定,Jacobi迭代方法收敛.

由（1）知，当时，对称正定当且仅当时正定同理可得，所以，当时Jacobi迭代方法收敛，否则发散.96

松弛法SORSuccessiveoverrelaxation换个角度看Gauss-Seidel方法：其中ri(k+1)=相当于在的基础上加个余项生成.下面令，希望通过选取合适的来加速收敛，这就是松弛法.iikikikiarxx)1()()1(+++=w0<<1低松弛法

=1Gauss-Seidel法>1(渐次)超松弛法SOR

97

写成矩阵形式：松弛迭代阵定理

设A

可逆，且aii0，松弛法从任意出发对某个收敛

(H)<1。98

定理

（Kahan必要条件）设A

可逆，且aii0，松弛法从任意出发收敛0<<2.证明：

利用，收敛|i|<1总成立|det(H)|<1=|det(H)|=|1|n

<10<<2.

99Thm若系数阵严格对角占优，若0<

≤1

，则松弛法收敛.(P171)100

定理

（Ostrowski-Reich充分条件）若A对称正定，且有0<<2，则松弛法从任意出发收敛。(参阅冯康等编《数值计算方法》)Q:什么因素决定收敛的速度?考察迭代：设B有特征根1、…、n

对应n个线性无关的特征向量。则从任意出发，可表为的线性组合，即~Ans:

(B)越小，迭代收敛越快.对于SOR法，希望找使得

(H)

最小.101

定理

若A

为对称正定三对角阵，SOR的最佳松弛因子(opt)为

此时.例：，考虑迭代格式问：取何值可使迭代收敛？

取何值时迭代收敛最快？解：考察B=I+A

的特征根1=1+,2=1+3

收敛要求

(B)<1-2/3<<0(B)=

max{|1+

|,|1+3

|}

当取何值时最小？-2/3-1/30=-1/2102SOR法化为G-S迭代法G-S法为SOR法的特例,SOR法为G-S法的加速例用G-S法和SOR法求下列方程组的解,要求精度1e-6一个例子103Solve:(1)G-S迭代法104

1110.75000000.37500001.50000000.56250000.53125001.54166670.65104170.59635421.61458330.70182290.65820311.6727431……….0.99999330.99999231.99999260.99999430.99999351.99999370.99999520.99999441.9999946k=71x=0.9999950.9999941.999995满足精度的解迭代次数为71次105(2)SOR迭代法1110.63750000.01218751.31990630.20042700.37175721.31228050.65503350.53401191.69228480.70584680.77334011.7771932………..0.99999900.99999761.99999910.99999840.99999931.99999890.99999980.99999941.999999