数值计算方法第三章

第三章线性代数方程组的解法针对有惟一解的非齐次线性代数方程组求解问题，本章主要主要介绍Gauss(高斯)顺序消去法、主元素消去法、三角分解法等直接求解算法，以及Jacobi(雅可比)法、Gauss-Seidel(赛德尔)法和逐次超松弛等迭代求解算法.§3.1弓I言大量的科学与工程实际问题常常可以归结为求解含有多个未知量土，码，…，土的线性代数方程组ax+axHFax=b1111221nn1ax+axFFax=b<2112222nn2(3.1)………ax+axFFax=bn11n22nnnn其对应的矩阵形式为Ax=b，其中n阶非奇异矩阵A以及n维列向量x和b分别定义如下A=(a11a21a12a22…a1n'a2nx=、x2b=rb)人1b2a•••aa:\x\:,b,V・n1n2.•..nn)nVnJ线性代数方程组数值解法可以分为直接法和迭代法两类.所谓直接法，是指在没有舍入误差的假设下，经过有限步运算就能得到方程组精确解的一类方法；而迭代法则是从一个初始向量x(0)出发，按照一定的迭代格式产生一个向量序列{X(k)}+8，当该序列收敛时，其极限就是线性方程组的解.k=0§3.2Gauss消去法Gauss消去法的基本思想是通过逐步消元将线性方程组(3.1)化为系数矩阵为三角形矩阵的同解方程组，然后用回代算法解此三角形方程组，并得到原方程组的解.一、三角形方程组解法对于如下形式的下三角形方程组，

TOC\o"1-5"\h\z"1=b1(3.2)<ax+ax-b2112222ax+axHbax—bn11n22nnnn若a.^0(i-1,2,,n)，则(3.2)的解为(3.2)ll<x1-b\an7(3.3)x—(b-ax-ax--ax)/a(k—2,3,...,n)Ikkk11k22k,k-1k-1kk此方法称为前推算法・...类似地，对于如下形式的上三角形方程组，ax+ax+111122aax+ax+111122a22x2++ax—b1-b21nn+ax(3.4)若a..壬0(i-1,2,,n)，贝U(3.4)的解为x-乂

nann〔、k=ka,k+1xk+1Un)/akk(k-n-1,n-2,,1)(3.5)此方法称为回代算法・……由于求解三角形方程组的过程很简单，因此只要把方程组化为等价的三角形方程组，求解过程就很容易完成.二、Gauss顺序消去法设线性方程组（3.1）的系数矩阵A-气）nxn非奇异，为描述方便，记该线性方程组的增广矩阵为m)=(、程组的增广矩阵为m)=(、a(1)a(1)…a⑴a⑴11121n1,n+1a⑴a⑴…a⑴a⑴21222n2,n+1其中……1a⑴a⑴…a⑴a⑴1'n1n2nnn,n+1/(3.6)a(i)-a,a(i)-b(i,j-1,2,...,n)(3.7)Gauss顺序消去法实际上就是在求解过程中方程次序不变（即没有换行操作）的Gauss消去法，具体过程如下：

TOC\o"1-5"\h\z第1步设a(1)丰0，用-“乘以(3.6)第1行后加到第，行(，二2,3,,n),\o"CurrentDocument"11a⑴11则…第，•行的第j个元素化为此时增广矩1a(1)a(1)-^^-a(1)=:a(ija⑴ijj11阵(3.6)相应地化为a(1)a(1)…a⑴a⑴11121n1,n+1a(2)…a(2)a(2)222n2,n+1…2)(j=2,3,n十•••1(3.9)(3.8)ka(2)…n2a(2)nna(2)n,n+1/TOC\o"1-5"\h\z第2步设a(2)莉，用-—乘以(3.9)第2行后加到第，行(i=3,4,...,n),22a(2)22则第，•行的第j个元素化为(3.10)a⑵一a-a⑵=:a(3)(j=3,4,n十122增广矩阵(3.9)相应地化为…ra(1)a(1)(3.10)ra(1)a(1)a⑴…a(1)a(1)1112131n1,n+1a(2)a(2)…a⑵a(2)22232n2,n+1a(3)…a(3)a(3)333n3,n+1:•••...•.••(a(3)a(3)a⑶1n3nnn,n+1(3.11)(3.12)类似地，当完成了第1〜n-1步消元后，(3.11)就化为(3.12)ra⑴a⑴…a⑴a(1)11121n1,n+1a(2)…a(2)a(2)222n.2,n+1.•…Ia(n)a(n)nnn,n+1/此时已把原方程组(3.1)等价转化成了上三角方程组(3.12)，可用回代算法求解该上三角方程组，回代公式为a(n)x—n,n+1na(n)na(n)nn1x—ka(k)kk(3.13)na(k)-乙a(k)xk,n+1kjjj—k+1(k—n-1,n-2,,1)由上面的消元过程知,Gauss顺序消去法的消元过程能进行下去的前提条件是akk)丰0(k—1,2,,n-1).回代过程则要求akk)丰0(k—1,2,,n).例3.1用Gauss顺序消去法解线性方程组（计算过程保留4位小数）0.0003x+6x—4.0001<126x-9x—-412解将第一个方程乘以-6再除以0.0003然后加到第二个方程得'0.0003x+6.0000x—4.0001<-120009.0000：—-80006.00002可得x2—0.6667,代入方程第一个方程得x1—0,而方程组的精确解为x1—3，x—2，可见用Gauss顺序消去法得到的解完全失真.23例3.1说明，当a（k）壬0（k—1,2,,n）时，Gauss顺序消去法能够进行下去，但当|。修）|-0或时？|相对于其他元素很小时，消去过程在浮点计算过程中产生的舍入误差将导致计算结果的误差急剧增大・在这种情况下，可以采用Gauss主元消去法.三、Gauss主元消去法对于上面的例3.1，如果我们先交换两个方程的顺序，得到等价方程组6x-9x—-4〔0.0006x1+6x2—4.0002计算过程中同样保留4位小数，用Gauss顺序消去法可解得x2—0.6667，气—0.3333.

上面的求解过程实际上就是Gauss主元消去法.根据主元选取范围的不同，Gauss主元消去法又分为列主元消去法和全主元消去法.1Gauss列主元消去法Gauss列主元消去法的执行过程如下：首先’在增广矩阵(3.6)的第1列元素中选绝对值最大的元素%,称之为第1列的主元，即有a⑴=maxa⑴Vgi1如果匕主1,则交换增广矩阵中第1行和第*行的元素，交换后的增广矩阵仍记为式(3.6)，但此时aa)已是第1列的主元.用主元aa)将其下边的n-1个元素a⑴(i=2,3,...,n)消元为零的过程与Gauss顺序消去过程的第一步相同，从而得增i1广矩阵(3.9).接着，在式(3.9)的第2列中除第1个元素外的其余n-1个元素中选主元am,即有a（2）i22=maxa（2）2<i<ni2如果iq,则交换(3.9)中第2行和第i行的元素，此时新的aa（2）i22=maxa（2）2<i<ni2当完成第1~n-1步的选列主元及相应高斯消去过程后，则得增广矩阵(3.12)，最后利用回代公式(3.13)求得原方程组(3.1)的解.列主元消去法除了每步需要按列选主元并作相应的行交换外，其消去过程与Gauss顺序消去法的消去过程相同.归纳起来有如下算法.算法3.1(列主元消去算法)输入方程组阶数小二维数组A存放的增广矩阵(AIb)；输出方程组的解Xj(i=1,2,,n)；Step1选主元的消去过程…对k=1,2,…,n-1

找行号+1,.../},使橹maxa(k)；k<i<n*若ik必则交换增广矩阵A的第k行和第七行元素;消元：对i=k+1,k找行号+1,.../},使橹maxa(k)；k<i<n*ikkjij又对j=k+1,k+2,•••,n+1，计a(k+D=a(k)—la(k)；ijStep2回代过程a(n)—n,n+1a(n)nn(2)对k=n-1,n-2,ikkjij(2)对k=n-1,n-2,…,1，a(k)

kka(k)—和a(k)xj=k+1不同于Gauss顺序消去法，Gauss列主元消去法只要线性方程组的系数矩阵A非奇异，即可求出其惟一解.2Gauss全主元消去法全主元消去法选主元的范围更大.对增广矩阵(3.6)来说，第1步是在整个系数矩阵中选主元，即选绝对值最大的元素，经过行列交换将其放在。口元素的位置，然后实施第一步消元过程.第二步是在(3.9)中的n—1阶矩阵a(2)…a(2)・22….2na(2)a(2)n2nn中选主元，其余各步选主元的范围依次类推，而每步选主元后的消去过程同列主元法的消去过程.全主元消去法每步所选主元的绝对值，通常比列主元消去法所选主元的绝对值要大，因而全主元消去法的求解结果更加可靠.但全主元消去法每步选主元要花费更多的机时，并且由于对增广矩阵可能进行列交换，使得未知量X,x，…,x的12n次序发生了变化，在编程时应当注意.§3.3矩阵三角分解法前面3.2节中介绍过，对于下三角形方程组和上三角形方程组很容易通过前推算法和回代算法求解.如果能将线性方程组Ax=b的系数矩阵A分解成A=LU(其中L是下三角阵，U是上三角阵)，则通过求解下三角方程组Ly=b得到列向量y，再通过求解上三角方程组Ux=y即可得到原方程组的解x・一、矩阵三角分解的基本概念Gauss顺序消去过程实际上就是对增广矩阵(AIb)进行若干次初等行变换，使之化为(UIg)的形式，其中U为上三角矩阵.由线性代数理论知，对一个矩阵进行一次初等行变换，相当于给这个矩阵左乘一个相应的初等矩阵.定义如下的一组单位下三角矩阵(对角线为1的下三角矩阵)1-1001Ln1-一1・.•11Lk=-1k+1,k-101n,k其中...010—1321•••<•••0—120「1-1Ln-1=1—1n,n—11=ik(k=1,2,.n—,ika(k)kk则对应的消元过程依次为i=;k+k4,2n,Ln-1L1[A(i)0(1)]=[A⑵|b⑵L「A(k)|b(k)]=「A(k+1)|b(k+1)A(n-1)|b(n—1)]=「A(n)|b(n)]=：[U|g]归纳消去过程有L1L2L2L1(A|b)=(U|g)(3.14)由(3.14)得...L1L2L2L1A=U(3.15)记L=(LLLL)-1(L仍是单位下三角阵)，则有n-1n-221…A=LU,Lg=b(3.16)定义3.1若n(n>2)阶矩阵A可以分解为一个下三角阵L和一个上三角阵U的乘积，即A=LU，则称这种分解为矩阵A的一个三角分解或LU分解；若L是单位下三角阵，U是上三角阵，则称为Doolittle(杜里特尔)分解；如果L是下三角阵，U是单位上三角阵，则称为Crout(克劳特)分解.定理3.1[26]如果n阶(n>2)矩阵A的前n-1阶顺序主子式均不为零，则A有惟一Doolittle分解和惟一Crout分解.下面我们以Doolittle分解为例说明矩阵三角分解的具体过程.二、Doolittle分解设矩阵A的所有顺序主子式不为零，并且有Doolittle分解A=LU，则求解线性方程组Ax=b等价于求解下三角方程组Ly=b及上三角方程组Ux=y,即'1l21(b1]

b2(3.17)"'nl'n2uu1112'1l21(b1]

b2(3.17)"'nl'n2uu1112u22u1nu2n用前推算法和回代算法，易求得这两个三角形方程组的解分别为'1=气'=b-Ul'(k=2,3,,n)kkkjjj=1(3.18)(3.19)(3.20),1)x=—(y-ux)(k=n-1,n-2,kukkjj1kkj=k+1从而得到原方程组的解x=(气,x2，...,X〃)T.a11a12a11a12.••a1n'\1(u11u12.••u)1na21•••a22•,••••a2n•,•=l21.••1...•••u22•••u2n•••iaVn1an2.••a/nniln1ln2.••1JVuJnn下面我们从矩阵乘法出发，推导Doolittle分解的具体计算公式.设A=LU,即第1步根据矩阵乘法规则，对A的第一行兀素和第一列兀素有(3.21)(3.22)(3.23)-Eikjm=1ukml=』(a-Elikuikkkm=1(j=kk+1,n,mj•••(3.22)(3.23)-Eikjm=1ukml=』(a-Elikuikkkm=1(j=kk+1,n,mj•••u)i(=k+k+immk)(k=2,3,…,n)2,n,)(3.25)u=a(j=1,2,,n)1j1jal=—i1(i=2,3,••;n)z'1u11一般地，设矩阵U的前k-1行和矩阵L的前k-1列的元素已经求出.第k步由矩阵乘法，让第k行元素、第k列元素对应相等，得a=Elu+u(=,k+1,)nkjkmmjkjSm=1(3.24)a=Elu+lu(=+,+2,)nikimmkikkkm=1于是求得矩阵U的第k行和矩阵L的第k列的元素为式(3.23)和式(3.25)就是矩阵A的Doolittle分解的计算公式，其计算特点是：U和L的元素一行一列交叉进行，先求矩阵U的行元素，再求L的列元素；如图3.1所示的逐框进行.

n1n2n1un2unn第n1n2n1un2unn第n框图3.1Doolittle分解计算次序如果每步将计算结果u和/仍然存放在矩阵A的相应元素。和。所占的单元kjikkjik上，不再占用新的单元，这种存贮形式称为紧凑格式.比较式式(3.29)和式(3.25)只需直接对增广距阵进行分的第一个公式’可以看出，七和u.j的计算方法完全一致’因此实际按紧凑格式计只需直接对增广距阵进行分aaaabuuuuy1112131411121勺aaaabluuuy212223242—^21—22-123242aaaablluuy3132333433132—33-134’3aaaabkllluy)414243444/I4142434414k例3.2用直接三角分解法解线性方程组算时，无须对线性方程组的右端向量进行单独处理，解即可.以四阶矩阵为例，可表示如下+3x=342x+x+2x=1x一x+2x+2x=32x+2x+5x+9x=71234112302121-1222259(10112302121-1222259(1000101-11k201解系数矩阵A=的LU分解为kL=/0)0017(100k012002110^3、/3、/3、/1、由Ly=1，解得y=1，再由Ux=1，解得x=0311117J0[UJ0[UJ10J求解过程用紧凑格式表示为[012°211-12、22533[112102T231-19j712021113132111；°回代公式(3.20)求得所给方程组的解0八--1八x==0,x=1一x=1,x=(1-2x一x)=0,x=3一3x一2x一x=1423422431432实际应用中对阶数较高的线性方程组也应该用选主元的三角分解法求解，以确保计算结果的可靠性，具体求解过程参见参考文献[6].三、平方根法许多实际问题所归结的线性方程组的系数矩阵往往是对称正定的，对于这种类型的线性方程组可建立计算量更小的三角分解法.定义3.2设A是n阶(n>2)对称正定矩阵.称A=LLt为矩阵A的Cholesky(乔列斯基)分解，其中L是非奇异的下三角阵.定理3.2[27]n阶(n>2)对称正定矩阵A一定有Cholesky分解A=LLt.当限定L的对角元全为正时，该分解是惟一的.下边推导Cholesky分解的计算公式.由A=LLt知[aa1112aa2122a[aa1112aa2122a1na2n[l11l21l22丫l11l21…ln11l22…L(3.26)\aan1n2annfn1'n2lnn第1步由矩阵乘法知a=l2,a=l-1，进而求得1111i1i111(3.27)l=寸a,l=la(i=2,3,,n)(3.27)ii

一般地，设矩阵L的前k-1列元素已经求出.第k步由矩阵乘法知a-区12a-区12+12,

kkkmkkm-1aik与’i"/m-1(i-k+1,k+2,,n)(3.28)于是疽’mm-疽’mm-1l--1(a-芸lliklikimkmkkm-1lkk(i-k+1,k+2,(k-2,3,,n),n)(3.29)式(3.27)和式(3.29)就是对称正定矩阵•A的Cholesky分解计算公式.利用Cholesky分解法求解方程组Ax-b时，由于分解公式中有开方运算，故也称Cholesky分解法为平方根法・例3.3用平方根法求解线性方程组(41-10)(x1:x(7)13-1082—-1-152x-43V0024,VxJv6/由Cholesky分解可得(41-1V01(41-1V013-10-1-1520)024121——202324!2<2251——232寸亘求解下三角方程组121121——20232*打('y2^3(7)-8-4Vyjv6J得3,y,y,y)=0225；2、石,-6..以丘,2、.商⑸t，再求解上三角方程组1234''''12血212血21——232<亘(气］

尤2尤3气)7)225

次62屈

kJ得(x,x,x,x)t=(1,2,-1,2)t.1234四、解三对角方程组的追赶法在二阶常微分方程边值问题、热传导问题以及三次样条插值等问题的求解时，经常要求解系数矩阵严格对角占优的三对角线性方程组Ax=d，即(bc11(bc11abc222・...abn—1n—1ancn—1bnd1

d•2:dkdn)(3.30)其系数矩阵A是严格对角占优的，即吟|>CJ(3.31)<\bI>a+c(i=2,3,,n—1)」b」>lanl...(3.31)根据线性代数理论，当矩阵A严格对角占优时，其各阶顺序主子式必不为零，故线性方程组(3.30)的系数矩阵A必有惟一的Doolittle分解.由于矩阵A的三对角特点，矩阵A有如下更特殊的三角分解形式(bc)(111(bc)(111abcl2222・・・・•••abcn—1n—1n—1kab)knn1l1n—1ln、(uuc1u2c2・・un—1cn—1u)1八按照矩阵乘法由式(3.32)可得如下三角分解公式L'<匕=兽(i=2,3,,〃)(3.33)，Tu=b-1c…1iiii-1这样，求解方程组Ax=d就转化为求解两个三对角方程组Ly=d和Ux=y,其计算公式为]'广《(3.34)[y=d-ly'(i=2,3,,n)nu5x=一(y一cx)(=n-1,n-2.,1)iUii+1iIi(3.35)计算过程(3.33)〜(3.35)称为解三对角方程组(3.32)的追赶法或Thomas(托马斯)方法.§3.4解线性方程组的迭代法本节介绍求解线性代数方程组的迭代法.首先给出向量范数和矩阵范数的概念.一、向量与矩阵的范数定义3.3设定义在n维向量空间Rn上的非负实值函数||•||满足正定性：IIXI20,VxGRn.当且仅当x=0时，||x|=0;齐次性：忡x||=|a|||x||,VxGRn,aGR;三角不等式：||x+y||<|x||+||y||,Vx,ygRn；则称实值函数||・||为向量空间Rn上的一种向量范数・对任意向量x=(x1,x2,,xn)TGRn,分别定义如下实值函数

£政1i=12i=1=maxx1<i<n1可以验证，它们分别都满足定义3.3的三个条件，因而都是向量范数，这是常用的三种向量范数，分别称为2i=1=maxx1<i<n1定义3.4设｛x俄）声是Rn中的向量序列，若有向量xyRn，使k=0limg）-x*||=0，则称向量序列｛x（k）｝+8收敛于向量x*，记为limx«）=x*.k—8k=0k—8（1）正定性：||A||>0,VAeRnxn.当且仅当A=0时，||A||=0;（2）齐次性：|aA||（1）正定性：||A||>0,VAeRnxn.当且仅当A=0时，||A||=0;（2）齐次性：|aA||=|a|||A|，VAeRnxn,aeR；（3）三角不等式：||A+B||<||A||+||B||,VA,BeRnxn；（4）自相容性：||AB||<||A||||B||，VA,BeRnxn.则称||・||为Rnxn上的矩阵范数.设A=、），常用的矩阵范数有||A||=max^|（列范数）1<j<f"||A||=3（A「A）（谱范数）、maxIIAIIIIAIIF=max（行范数）ij切1<i<ni=12ij（F范数）定义3.6设nxn阶矩阵B在复数范围内的诸特征值为人（i=1,2,,n），则称ip（B）=max人|为矩阵B的谱半径.1<i<n'定理3・3对任意矩阵BeRnxn有p（B）<||B||，其中||・||是Rnxn上的任何一种矩阵范数.证明设人是B的任一特征值，x=0是相应的特征向量，则Bx=入x显然有向量yeRn，使xyT为非零矩阵.用yT右乘上式得

Bxyr=人xyT由矩阵范数定义，WxyvIBII-||xyBxyr=人xyT由矩阵范数定义，WxyvIBII-||xyt||，由此得^<||B||设有〃阶线性方程组(3.36)其中系数矩阵A非奇异，向量b^0，方程组(3.36)有惟一解x*.将方程组(3.36)等价变形为(3.37)(3.38)x=Bx+g并定义迭代过程如i)=Bx(k)+g(k=0,1,)k=0当取定初始向量x(0)GRn后，式(3.39)便产生一个向量序列{x(k)}+8,若它收敛于某向量x*，则x*一定是(3.37)的解，当然也是原方程组的(3.36)解•称形如(3.38)的迭代法为简单迭代法，并称B为该简单迭代法的迭代矩阵.(3.37)(3.38)k=0显然，方程组Ax=b的等价形式(3.38)不惟一，因而可建立不同的简单迭代法.，对相同的初始向量x(0)，产生的向量序列就会不同，有的可能收敛，有的可能不收敛.当收敛时，只要k充分大，就可以将x(k+i)作为方程组的近似解.另外，对同一简单迭代法(3.38)可能关于某个初始向量x(0)产生的序列收敛，而关于另外一个初始向量产生的序列又不收敛.在本节讨论中，只有关于任意初始向量x(0)都收敛时，我们才称迭代法收敛.定理3.4简单迭代法x(k+i)=Bx(k)+g(k=0,1,2,)对任意初始向量x(0)都收敛的充要条件是limBk=0....证明设x*是方程组的解向量，则有

(3.39)另外，根据迭代格式，有如=Bx(k-1)+g(k=1,2,)(3.40)由式(3.40)和式(3.39)得...(3.39)x(k)—x*=B(x(k-1)—x*)=B2(x(k-1)—x*)==Bk(x(0)-x*)(3.41)由式(3.41)可得，对任意初始向量x(0)迭代法收敛，即limx(k)=x*的充要条ks件是limBk=0.ks定理3.5[18]简单迭代法x(k+1)=Bx(k)+g(k=0,1,2,)对任意初始向量x(0)都收敛的充要条件是迭代矩阵的谱半径p(B)<1.…结合定理3.3及定理3.5，我们可得关于收敛性判定的如下充分条件.定理3.6设简单迭代法x妇)=Bx(k)+g(k=0,1,2,)的迭代矩阵B的某种矩阵范数小于1,即有||B||<1,则该简单迭代法关于任意初始向量x(0)收敛.矩阵范数|B|[和|B||容易计算，故实践中常用条件||B||]<1或||B||<1来判定简单迭代法的收敛性.下面介绍三种最常用的简单迭代法：Jacobi(雅克比)迭代法、Gauss-Seidel(高斯-赛德尔)迭代法以及逐次超松弛迭代法.三、Jacobi迭代法设n阶线性方程组(3.1)的系数矩阵A非奇异，且a^丰0(/=1,2,..，n),将方程组(3.1)改写成x=(b-ax-ax1122133x=(b-ax-ax2211233-axJ.a1nn':22(3.42)进而可写出简单迭代法(b—ax—ax—nn11n22-axn,n-1n-1nn

x(k+1=\b—ax炊—ax.k—)—axn1x(k+1=b—axk(tax=(b-ax-ax1122133x=(b-ax-ax2211233-axJ.a1nn':22(3.42)进而可写出简单迭代法(b—ax—ax—nn11n22-axn,n-1n-1nnx(k+1=\b—ax炊—ax.k—)—axn1(3.43)xk—)22Xkn—,n—式中a-^2aiia-^1a22a—^n1[aa■n2a%1a11aa22.0.=—D-(L+U)，gj=[，111ba22=D-ib(3.44)/01a'0a1112a11nL=a：210••D=a12•U=0•a--2n.:.n1•an2•0/L•annk•1.••0J7根据定理3.5和定理3.6,Jacobi迭代法关于任意初始向量x（0）都收敛的充要条件knnnnnn定理3.7设系数矩阵A=（a）严格对角占优，即a』斗ijnxn(i=1,2,,n)（按行严格对角占优）(3.45)或者"IJj=1J技•••laj>£Ki=1i丰J(j=1,2,,n)•••（按列严格对角占优）(3.46)成立，则求解Ax=b的Jacobi迭代法关于任意初始向量x（0）收敛.证明B==—D-1(L+U)=—D-1(A—D)

设人为B的某个特征值，则有det(人I—B)=det(人I+D-i(A—D))=det(D-i)det(人D+A—D)=0由detD-1丰0知(3.47)(3.48)(3.49)人D+A—D=fXa11a21a12Xa22••-:...a1n'a2n:若|X(3.47)(3.48)(3.49)人D+A—D=fXa11a21a12Xa22••-:...a1n'a2n:若|X>1,则由(3.45)可得•ka'n1••••an2♦♦♦-•XaJnnlxai,或由(3.46)可得>ai*(i1ijj=1j归(i=1,2,n,)•••而det(人D+A—D)=0(j=1,2,,n)|人七|>\aa〃」ea.i=1i=j关于迭代法的误差控制，设£为允许的绝对误差限，可以通过检验|X(k+1)—x(k』=maxx(k+1)—x(k)v£是否成立，来决定计算是否终止迭代・'"“算法3.2(Jacobi迭代法算法)输入方程组系数矩阵A，维数n，右端向量b，初始向量x(0)，精度要求£，最大迭代次数N;输出方程组的解向量x;Step1令k=1；Step2对i=1,2,,n依次计算1X⑴=——ia.ax(0)ax(0)ijjj=1j丰iStep3如果maxx（i）-x（o）<e，则输出x（1），算法结束；否则转Step4；Step4若k<N，则令k=k+1,x（0）=x（i），转Step2；否则输出迭代失败信息，算法结束.例3.4用Jacobi迭代法求如下方程组的近似解x（s），要求先讨论收敛性，若收敛，取x（0）=（0,0,0>,迭代至||x（k+i）-x（k）||<10-3.5x一x+x=10<x-10x一2x=27一x+2x+10x=13V123解由于系数矩阵按行严格对角占优，因此Jacobi迭代法收敛.构造Jacobi迭代格式为x（k+1）=0.x(k)-0x2+2123x（k+1）=0.x(k)—0x2一2.7213x（k+1）=0.x(k)-0x2)+1.3312取x（。）=（0,0,0>，计算结果如表3.1所示.表3.1Jacobi迭代法计算结果kx(k)1x(k)2x(k)312.0000000-2.70000001.300000021.2000000-2.76000002.040000031.0400000-2.98800001.972000041.0080000-2.99040002.001600051.0016000-2.99952001.998880061.0003200-2.99961602.000064071.0000640-2.99998081.9999552由于maxx⑺一x(6)<10-3，故取x*-x⑺=(1.0001,-3.0000,2.0000》.四、Gauss-Seidel迭代法及JGS迭代法如果简单迭代法(3.38)收敛，土(奸1)应该比x(k)更接近于原方程组的解x*(i=1,2,n)并且在计算第i个分量x(k+1)时，前面的i-1个分量已经算出，因ii而可以用新值x(k+1)，…，x(：1)来代替旧值x(k)，…，x*，希望能提高收敛速度.这样得到如下形式的迭代法xU+1)=旧bxQ+1)+Ebx(k)+g(i=1,2,,n)(3.50)j=1j=i称该迭代法为与简单迭代法(3.38)对应的Gauss-Seidel迭代法・若将简单迭代法(3.38)的迭代矩阵B=(b)^分解为B=B1+B2，其中(0)(bbb\11121nb0bbB=21B=22•--2n12b•n2…..0>bj•-nn则简单迭代法(3.38)可写成...x(k+1)=Bx(k)+Bx(k)+g(k=0,1,2,)(3.51)TOC\o"1-5"\h\z\o"CurrentDocument"12对应的Gauss-Seidel迭代法可表示为...x(k+1)=Bx(k+i)+Bx(k)+g(k=0,1,2,)(3.52)式(3.52)等价于如下简单迭代法...x(k+1)=(I-B)-1Bx(k)+(l-B)-1g(3.53)Gauss-Seidel迭代法的迭代矩阵为B®=(l-B,1B2.这样就可以使用简单迭代法收敛性的各种判别方法来判别Gauss-Seidel迭代法的收敛性.下面再给出Gauss-Seidel迭代法的一个收敛性判别定理.定理3.8设简单迭代法(3.51)的迭代矩阵B=B1+B2满足|B|<1或|B1<1，则相应的Gauss-Seidel迭代法(3.52)关于任意初始向量都收敛.证明仅以|B|=maXn|b」<1为例进行证明，将式(3.51)的第i个方程与yj=1x=Bx+g相应的方程相减，有

X(k+1)-X=(i=1,2,,n)旧bC(k+1)-x)+£bCx(k)—x)jTj=i记8=max成)一x,k1<i<nX(k+1)—X.设x(k+1)一xiX(k+1)-X=(i=1,2,,n)记8=max成)一x,k1<i<nX(k+1)—X.设x(k+1)一xi0于是有进而有a=第，C斗.1,则由上式得<S|bjj=1j=1j=x(k+1)—x,+,+£|b|x(k)—xj<8a+8P-jj=ik+1iki=8，由于式(3.54)对i=1,2,,n都成立，取i=i°得8k1<8a+8Pk+1i0‘0(3.54)8k+1<PE—T8k’0(3.55)0<a+P.=?勺<l|B|土<11—aiL=maxi(i=1,2,,n)(i=1,2,,n)(3.56)由式(3.55)和式(3.56)知8<L8<<L+180(3.57)故8k+1maxx(k+1)—xt0(kT3),艮即Gauss-Seidel迭代法(3.52)关于任意初1<i<n结合定理3.6及3.8知，当简单迭代法的迭代矩阵范数||B||<1或|B|<1时，简单迭代法(3.38)及与之相应的Gauss-Seidel迭代法(3.52)同时关于任意初始向量收敛.

卜面讨论与Jacobi迭代法相对应的Gauss-Seidel迭代法.可以将式(3.43)改写为x(k+1)1x(k+1)，2写为x(k+1)1x(k+1)，2—ax(k+i)—ax(k)—.—a__x(k+1)n,n—1n—1x(k+1)=\b一ax(k+1)一ax(k+i)一、,方••••"••••,分1—a__x(k+1)n,n—1n—1(3.58)1x(k+1)=一

「aViivlb—Eax(k+1)—乙x(k)'"尸i+1〃'/(i=1,2,・・,n)该迭代格式称为与Jacobi迭代法相对应的Gauss-Seidel迭代法，(3.59)简称JGS迭代法.JGS迭代法是一种特殊的Gauss-Seidel迭代法，因此关于Gauss-Seidel迭代法的收敛性判定定理均适用于JGS迭代法.若设系数矩阵A=L+D+U，其中L,D,U的含义同(3.44)式，可得JGS迭代法的迭代矩阵为(3.60)BJGS=-(D+L)-1U(3.60)定理3.9设系数矩阵A")nxn严格对角占优，则求解Ax=b的JGS迭代法对任意初始向量x(0)都收敛.类似于定理3.7的证明过程可以证得该定理.定理3.7和定理3.9说明，线性方程组的系数矩阵A严格对角占优时,Jacobi迭代法和JGS迭代法对任意初始向量都收敛.定理3.10[18]若线性方程组的系数矩阵A对称正定，则求解Ax=b的JGS迭代法对任意初始向量都收敛.算法3.3(Guass-Sediel迭代法算法)输入方程组系数矩阵A，维数n，右端向量b，初始向量x(0)，精度要求e，最大迭代次数N;输出方程组的解向量x;Step1令k=1；Step2对i=1,2,,n依次计算•1x⑴=iaiiVb—Eax(1)—Eax(0)(i=1,2,,n)*八j=i+1jj•1x⑴=iaiiVStep3如果max|x⑴-x(o)<s，则输出x⑴，算法结束；否则转Step4；Step4若k<N,k=k+1Step3如果max|x⑴-x(o)<s，则输出x⑴，算法结束；否则转Step4；例3・5用JGS迭代法求解例3.4中的线性方程组，取x(0)=(0,0,0》,迭代终止条件为|x(k+1)-X(k)[<10-3.8解由于系数矩阵按行严格对角占优，因此知JGS迭代法收敛.JGS迭代格式为x(k+1)=10.2(k)-20x2+32x(k+1)=0.X(k+1)—0X2(T2.7213x(k+1)=0.X(k+1—0x2+1)+1.31321取x(0)=(0,0,0》，迭代结果如表3.2所示.表3.2JGS迭代法计算结果kX(k)1X(k)2X(k)312.0000000-2.50000002.000000021.1000000-2.99000002.008000031.0004000-3.00156002.000352040.9996176-3.00010861.999983550.9999816-2.99999851.9999979由于已有max|x(5)-x(4)<10-3，故可取x(5)作为原方程组的解.一般情况下:JGS迭代法比Jacobi迭代法收敛快.但应该注意，并不是任何时侯JGS迭代法都比Jacobi迭代法收敛快，甚至有Jacobi迭代法收敛而JGS迭代法不收敛的例子.五、逐次超松弛迭代法设线性方程组Ax=b的系数矩阵A非奇异，且%主0(i=1,2,/),设已求得方程组的第k次近似x(k)，称r(k)=b-Ax(k)为第k次迭代的残差向量.当x(k)不

满足误差要求时，需要求第k+1次近似解xl+1).将残差向量r(k)乘以修正因子«D-1，补偿到x(k)，建立如下迭代x(kx(k+1)=x(k)+WD-1r(k)(k=0,1,2,…)(3.61)其分量形式为WX(kWX(k+1)=X(k)+—

//aiijXj)-f)Ij=1j=i)(3.62)(i=1,2,…,n;k=0,1,2,…)式中W称为松弛因子.考虑到计算X(k+1)时，X(k+1)(1<j<i-1)均已算出，于是将式TOC\o"1-5"\h\z(3.62)改写为'jX&+1)=X*)+仪a(3.63)b-国ax(k+1)-X&+1)=X*)+仪a(3.63)iijjijjj=1j=i(i=1,2,...,n;k=0,1,…)称(3.63)为逐次超松弛迭代法，简称SOR(SuccessiveOver-Relaxation)方法.式(3.63)还可整理为1X(k+1)=(11X(k+1)=(1-W)X(k)+Wiiaii\j"E)j=i+1)(b—Wax(k+1)-£】j=1(3.64)(i=1,2,—,n)当松弛因子w=1时，SOR即为JGS迭代法.因此，SOR方法也可看做JGS方法的计算值与当前近似向量球)的一种加权平均.如果松弛因子选取得适当，贝0SOR方法往往具有加速收敛的作用.若将系数矩阵A仍分解为A=L+D+U，则SOR方法的矩阵形式为x(k+1)=G-W)x(k)+wD-1(-Lx(k+1)-Ux(k))(3.65)与之等价的简单迭代法为(3.66)x(k+1)=Bx(k)+w(D+wL)-1b(3.66)其中B①=(D+®L、i61-①b-®U】(3.67)是SOR方法的迭代矩阵.关于SOR方法收敛性，除了可以用充要条件p(B)<1,充分条件||B||<1来判定外，还有如下定理.°°定理3.11[6]设AGRn.n且。济0(i=1,2,,n),则求解Ax=b的SOR方法对任意初始向量x(0)都收敛的必要条件牛是0<°<2.定理3.12[6]设矩阵a对称正定且0〈①七方，则求解Ax=b的SOR方法对任意初始向量x(0)收敛.定理3.13[6]设矩阵A严格对角占优且0<°V1,则求解Ax=b的SOR方法对任意初始向量x(0)收敛.例3.6取°=1.3，初始向量x(0)=(0,0,0>，用SOR方法解线性方程组-4-45"X)-4-45"X)"7、1X=一32kX3>3k370)5J要求X(+1)—X()<10-4.解系数矩阵A=解系数矩阵A="5-4<00、-1的顺序主子式气=5，5J气=9，气=40，且A=At，故系数矩阵A对称正定，所以0<°<2时SOR方法收敛.SOR方法的计算格式为°「_X(k+1)=(1一°)X(k)+——7+4X(k)15L2-<x(k+1)=(1-°)x(k)+兰-3+4x(k+1)+x(k)213」x(k+1)=(1—°)x(k)+—3+x(k+1)335L2」取店1.3,迭代结果如表3.3所示：表3.3SOR方法计算结果kX(k)X(k)X(k)12311.82000001.11280001.069328022.43131201.69274980.89931655332.85106621.91110621.007092642.95223061.97883200.9923685::::83.00018522.00016100.9999