数值分析 第6章 解线性方程组的迭代法

数理学院SCHOOLOFMATHEMATICSANDPHYSICS6.1引言6.2Jacobi迭代6.3Gauss－Seidel迭代6.4超松弛迭代法(SOR)第6章解线性方程组的迭代法对方程组bAx做等价变换gGxxbMNxMxbNxMxbxNMbAx11)(如：令NMA，则则可以构造序列gxGxkk)()1(若*)(xxkbAxgxGx***同时：*)(**)()()1(xxGGxGxxxkkk*)()0(1xxGk0kG所以，序列收敛与初值的选取无关.(6.1)（6.2）§6.1引言定义6.1：（收敛矩阵）0kG定理：矩阵G为收敛矩阵，当且仅当G的谱半径11)(0GGk证明：由于,)]([)(kkGG当,0limkkG故有,0)(limkkG即：.1)(,0)]([limGGkk反之，当1)(G时，由定理可知，对任意,0存在,||||G使,1)(||||GG于是：.0||||lim||||limkkkkGG从而.0limkkG,)(GG若有某种范数1pG，则迭代收敛.由于计算较困难，通常可利用)(G‖G‖＜1是判断收敛的充分条件，即使‖G‖＜1对任何范数都不成立，迭代序列仍可能收敛.§6.2Jacobi迭代nnnnnnnbxaxabxaxa1111111)(1)(1)(111,112232312122211212111nnnnnnnnnnnnbxaxaaxbxaxaxaaxbxaxaax)(1)(1)(1)(11,)(11)1(2)(2)(323)(12122)1(21)(1)(21211)1(1nknnnknnnknknnkkkknnkkbxaxaaxbxaxaxaaxbxaxaax写成一般格式：)(11)(11)()1(nijkjijijkjijiiikixaxabaxJacobi迭代算法1、输入系数矩阵A和向量b，和误差控制eps2、x1={0,0,…..,0},x2={1,1,…..,1}//赋初值3、while(||A*x2-b||eps){x1=x2;for(i=0;in;i++){x2[i]=0;for(j=0;ji;j++){x2[i]+=A[i][j]*x1[j]}for(j=i+1;jn;j++){x2[i]+=A[i][j]*x1[j]}x2[i]=(b[i]-x2[i])/A[i][i]}}4、输出解x2•迭代矩阵ULDA记nnaaD1100001,121nnnaaaL0000,1112nnnaaaU易知，Jacobi迭代有bxULD)(bxULDx)(bDxULDx11)(bDgADIULDG111,)(•收敛条件迭代格式收敛的充要条件是ρ(G)1，充分条件是||G||1.对于Jacobi迭代，例外还有一些保证收敛的充分条件.定义若满足nnijRaA)(niaanijjijii,,2,1|,|||1严格不等式成立则称A为严格对角占优矩阵；满足nnijRaA)(如果niaanijjijii,,2,1|,|||1若上式至少有一个严格不等式成立则称A为弱对角占优矩阵.定义假定,若存在排列矩阵使其中则称A可约，否则，若不存在排列阵P,使上式成立，则称A不可约.2,nRAnnnnRP2212110AAAAPPT)1(,)()(2211nrRARArnrnrr(6.3)•例如方程组123110011012321xxx就是可约的，求解时可先解2322121xxxx得解为，再代入最后一个方程得121xx03x若取它就是(6.3)的形式.210110011,001010100APPPT则定理若为严格对角占优矩阵或不可约弱对角占优矩阵，则且A非奇异.nnijRaA)(),,,1(0niaii(1)A为严格对角占优矩阵;(2)A为不可约弱对角占优矩阵.定理：若A满足下列条件之一，则Jacobi迭代收敛。例1.给定方程组3612363311420238321321321xxxxxxxxx它的精确解可构造如下迭代法若写成式(6.1)的形式，则迭代矩阵G及g可表示为：,)1,2,3(*Tx)3636(121)334(111)2023(81)(3)(1)1(3)(3)(1)1(2)(3)(2)1(1kkkkkkkkkxxxxxxxxx若取,按式(6.2)迭代10次可得误差它表明迭代序列收敛.12361133820,04121111011441830gG,)0,0,0()0(TxTx)999803.0,999838.1,000032.3()10(.000197.0||||*)10(xx程序如下：function[x,k]=Jacobi(A,b,x0,N,emg)n=length(A);x1=zeros(n,1);x2=zeros(n,1);x1=x0;k=0;r=max(abs(b-A*x1));whileremgfori=1:nsum=0;forj=1:nifi~=jsum=sum+A(i,j)*x1(j);endendx2(i)=(b(i)-sum)/A(i,i);endr=max(abs(x2-x1));x1=x2;k=k+1;ifkNdisp('迭代失败，返回');return;endendx=x1;A=[8-32;411-1;6312]；b=[20;33;36];X0=[0;0;0];[x,k]=Jacobi(A,b,X0,100,0.001)在命令窗口输入：x=3.000281567961441.999911821895700.99974047656463k=9运行结果：)(1)(1)(1)1(11,)1(11)1(2)(2)(323)1(12122)1(21)(1)(21211)1(1nknnnknnnknknnkkkknnkkbxaxaaxbxaxaxaaxbxaxaax§6.3Gauss－Seidel迭代)(11)(11)1()1(nijkjijijkjijiiikixaxabax在Jacobi迭代中，使用最新计算出的分量值Gauss-Seidel迭代算法1、输入系数矩阵A和向量b，和误差控制eps2、x2={1,1,…..,1}//赋初值3、while(||A*x2-b||eps){for(i=0;in;i++){for(j=0;ji;j++){x2[i]+=A[i][j]*x2[j]}for(j=i+1;jn;j++){x2[i]+=A[i][j]*x2[j]}x2[i]=(b[i]-x2[i])/A[i][i]}}4、输出解x2•迭代矩阵)()()1(1)1(bUxLxDxkkkbDUxDxLDIkk1)(1)1(1)(bDLDIUxDLDIxkk111)(111)1()()(bLDUxLDxkk1)(1)1()()(bLDgULDG11)(,)(是否是原来的方程的解？A=(D-L)-U系数矩阵对应的分解公式：•收敛条件迭代格式收敛的充要条件是ρ(G)1，充分条件是||G||1.另外还有一些充分条件定理：若A满足下列条件之一，则Gauss-Seidel迭代收敛。①A为严格对角占优矩阵②A为不可约弱对角占优矩阵证明：设G的特征多项式为)(sP，则ULDLDULDIGIPs)()()()(11ULDA为对角占优阵，则1时ULD)(为对角占优阵0)(ULD即0)(sP1即1)(G＃证毕定理若(1)解方程组Ax=b的GS法收敛.(2)若2D-A也对称正定，则J法也收敛.nnijRaA)(对称正定，则注：二种方法都存在收敛性问题。有例子表明：Gauss-Seidel法收敛时，Jacobi法可能不收敛；而Jacobi法收敛时，Gauss-Seidel法也可能不收敛。例2用GS法求例1中方程组的解.解GS法的迭代公式为取迭代5次得到)3636(121)433(111)2320(81)1(2)1(1)1(3)(3)1(1)1(2)(3)(2)1(1kkkkkkkkkxxxxxxxxx,)0,0,0()0(Tx,)000061.1,0000721.2,999843.2()5(Tx000157.0||||||||*)5()5(xxe这个结果与用J法迭代10次的误差相当.在一定条件下，若J法与GS法均收敛，则GS法比J法约快一倍，但也可能J法收敛而GS法不收敛或相反.function[y,n]=seidel(a,b,x0,emg)D=diag(diag(a));U=-triu(a,1);L=-tril(a,-1);G=(D-L)\U;f=(D-L)\b;y=G*x0+f;n=1;whilenorm(y-x0)emgx0=y;y=G*x0+f;n=n+1;endynA=[8-32;411-1;6312]；b=[20;33;36];X0=[0;0;0];[x,k]=Jacobi(A,b,X0,0.001)在命令窗口输入：运行结果：y=2.999842386641112.000072133594331.00006077328086n=5(0)1807109,8,(0,0,0)'9117Abx1、预处理9117180,71098Ab2、迭代格式例3)8(91)7(81)7(91)1(1)1(3)1(1)1(2)(3)(2)1(1kkkkkkkxxxxxxx(1)(2)(3)(4)(0.7778,0.9722,0.9753)'(0.9942,0.9993,0.9994)'(0.9999,0.9999,0.9999)'(1.0000,1.0000,1.0000)'xxxx3、结果101/21/2()1011/21/20BDLU1、Jacobi迭代211111112A特征值为3504IB12,350,2i101/21/2()01/21/2001/2BDLU2、Gauss－Siedel迭代12,310,2例4例5方程组Ax=b中，证明当时GS法收敛，而J法只在时才收敛.111aaaaaaA12/1a2/12/1a解只要证12/1a时A正定即可，由顺序主子式得于是得到12/1a1||0111,01221aaaa0)21()1(321det2233aaaaA2/1a对J法，由于迭代矩阵是J法收敛的充要条件.故J法只在时才