研究线性方程组迭代收敛速度.

1、研究解线性方程组迭代收敛速度一实验目的科学研究与生产实践中许多问题都可归结为线性方程组的求解，高效求解线性方程组成为了许多科学与工程计算的核心迭代法就是用某种极限过程去逼近线性方程组精确解的方法，该方法具有对计算机的存贮单元需求少，程序计算简单，原始系数矩阵在计算过程中不变等优点，是求解大型稀疏矩阵方程组的重要方法。常用的迭代法有Jacobi迭代法、Gauss一seidel迭代法、逐次超松驰法（SOR法）等。二实验摘要由迭代法平均收敛速度与渐进收敛速度的关系引入近似估计法，即通过对迭代平均收敛速度取对数，然后利用Mathematica软件对其进行拟合，给出拟合函数，最终得到了Jacobi迭代法

2、、Gauss一seidel法的平均收敛速度收敛到渐进收敛速度的近似收敛阶，以及逐次超松驰法（SOR法）的渐进收敛速度，且该法适用于其他迭代法收敛速度的估计。三迭代法原理1.Jacobi迭代法（J法）设方程组Ax=b，其中，A=（a.）GR"Xn，x,bGRn。jnxnA为可逆矩阵，可分裂为A=LU,其中，a21a31a32an1an2a0n,n-1a120a13a23a1na2n0an-1,n其中，BG称为式(1)的GaussSeide1迭代法的迭代矩阵。a11D二a22ann从而由Ax=b得到,xD-1(L+U)x+D-ib=D-1(DA)x+D-1b(ID-1A)x+D-1b令B

3、ID1A,fD1b,JJ由此可构造出迭代公式：x(k+1)Bx(k)+fJJ令初始向量x(0,0,0)，即可得到迭代序列，从而逼近方程组的解这种方法称为Jacobi迭代法,其中BJ称为Jacobi迭代矩阵。2.Gauss-Seidel迭代法(GS法)与Jacobi迭代法类似，将方程组Axb中的系数矩阵A分裂为AL+D+U,，其中D,L,U与前面相同。与Jacobi迭代法所不同的是，Gauss-Seidel迭代法将Jacobi迭代公式中的Dx(k+1)Lx(k)Ux(k)+b改为Dx(k+1)Lx(k+1)Ux(k)+b从而Axb可写成矩阵形式(L+D)x(k+1)Ux(k)+b，若设(L+D)

4、-1存在，则x(k+1)(L+D)1Ux(k)+(L+D)1b，其中，B(L+D)1U,f(L+D)1b，G于是GaussSeide1迭代公式的矩阵形式为x(k+1)BGx(k)+f。注：由于Gauss-Seide1迭代充分利用了迭代过程的新信息，一般来说，它的迭代效果要比Jacobi迭代好。3逐次超松弛方法(SOR法)根据Gauss-Seidel迭代法的迭代原理x仏+1)=Dt(Lx仏+1)Ux(k)+b)，我们将其修改为x(k+i)=(1w)x(k)+wX""D,即对x(k)和x"'"加入一个权因子，在这里我们称w为迭代公式的松弛系数。令X)

5、=D-1(Lx(k+1)一Ux(k)+b)(Gauss一Seidel迭代法)，贝U(k+1)x(k+1)=(1w)x(k)+wx=(1w)x(k)+wD1(Lx(k+1)Ux(k)+b)从而x(k+1)=(D+wL)1(1w)DwUx(k)+w(D+wL)1bG=(D+wL)1(1w)DwU，f=w(D+wL)1ww所以x(k+1)=Gx(k)+fww将其写成分量的形式严=(1-亦严+巴©-号肘刊-乞叫普*"2用aa匸1m+i我们称该公式为基于Gauss一Seidel迭代下的超松弛迭代公式。注：1)关于松弛系数w的选取，我们有以下性质：(i) .设方程组Ax=b，其中，A=

6、(a.)&Rnxn，x，bGRn，则解方程的SORjnxn迭代法收敛的必要条件是o<w<2；(ii) 若A为正定的对称矩阵，且0<w<2，则方程组Ax=b关于SOR迭代法是收敛的；(iii) 若A为正定的对称三对角矩阵，B为Jacobi迭代矩阵，若竹的谱半径P(Bj)<1，记p=p(Bj)，则SOR法的最优松弛因子为2W=一bl+°k且wp+Jw2屮-4(w一1)2/40<w<wp(G)=Sybww-1w<w<2b2) 由松弛系数的性质可知，通常只有当方程组的系数矩阵A为正定的对称矩阵时我们才使用SOR法；3) 当松弛系数

7、w二1时，SOR法记为GS方法；4) 关于(iii)提到的谱半径定义为：假设A为n阶可逆矩阵，九,九，九是它的n个特征值，我们称12np(A)=maxI九I1<i<ni为A的谱半径。容易证明A的谱半径是A的所有诱导范数的下确界，即P(A)二infIIAIIII-II其中，矩阵的范数如：IIAII二max£Ia1,11AII二max£Ia1,11AII二jp(AtA)。11<j<ni=1ijs1<i<nj=i可2下面我们就来讨论以上方法的迭代收敛速度，首先我们介绍收敛速度快慢的比较方法。四迭代法收敛速度的比较1.这两种迭代方法收敛性与BkT

8、0(kTs)是否成立有关，且收敛速度与Bkt0的速度有关。当p(B)<1时，Bk趋于零矩阵的速度有赖于p(B)的大小。一般说来，p(B)越小，则Bk趋于零矩阵的速度越快，反之就越慢。通常，当p(B)<1时，可以用正数-Inp(B)的大小作为迭代法渐进收敛速度的度量。这时-Inp(B)越大，迭代法的收敛速度愈大。2.平均收敛速度与渐进收敛速度之间的联系对于收敛的迭代法x二Bx+f（k二0,1,2.）,R二ln（llBkIh/k）称为平均收敛k+1kk速度（它与所用的范数以及k有关）；R=inp（B）称为渐进收敛速度。可以证明8RR=iimR，因此当kT8时假设成立下列渐近关系式古TC

9、（常数），为估计8kT8kRpk平均收敛速度收敛到渐进收敛速度的阶，当k充分大时有如下近似形式：Rk11UC,k>>1Rk两边取对数得：InRupinR+Inckk此式说明当k比较大时，inR与inR有近似的线性关系，而它们的斜率刚好k+1k为收敛阶P,这样可以通过当k比较大时作出inR与inR的拟合曲线来估计出Pk+1k值。五实验举例例1：考虑五阶方程组5x一3x一x=2145x+4xx=3125<2xx=11.其34x+2x=14150000.60.20.250000.250000.500.250000.50000.50、BJx+4x2x=0145渐进收敛速度为：R=ln

10、(p(B)=ln(0.66165261)=0.413014628J则迭代平均收敛速度R（见表1）以及取对数后作最小二乘拟合图像（见图1）如k下所示：表1kRklnRK10.22314355-1.4999400020.29891850-1.2075843030.35473695-1.0363787040.35677909-1.0306385050.37710417-0.9752338260.37622879-0.9775578370.38679290-0.9498658680.38586107-0.95227789图1即ln(R)=0.44226443ln(R)-0.53345499k+1k从而

11、得到收敛阶为p=0.442264432.该方程组的GaussSeidel迭代矩阵为：<0000.60.2、0000.150.3B二0000.50G0000.150.5510000.0750.275丿其渐进收敛速度为R=-ln(p(B)=-ln(0.425)=0.85566611gG则迭代平均收敛速度R（见表2）以及取对数后作最小二乘拟合图像（见图2）如下所k示：表2kRklnRK10.22314355-1.4999400020.35667494-1.0309304030.52300533-0.6481636240.60617053-0.5005939350.65606964-0.4214

12、883360.68933572-0.3720268770.71309721-0.3381375380.73091832-0.31345356图2ln（R）=0.584721431n（R）-0.11593362k+1k从而得到收敛阶为p=0.58472143。注：在本例中，由于方程组的系数矩阵严格对角占优，故前述两种迭代过程均收敛依实际迭代过程：对Jacobi迭代有"X（32）-X二6.4345694x10-7<10-6|x32|而GaussSeidel迭代有"X"-X""二5.5402987x10-7<10-6，即二者均收敛,|x1

13、8|但后者更快一些。这与收敛阶p=0.585>0.442之间关系一致。例2：考虑方程组4x+3x二2412v3x+4x一x二30123一x+4x二-2423用SOR迭代公式可得彳严=（1-町谭+鬟0-张严+铲）严=（1-亦缪+中-24+V）取初试量为x（o）=（0,0,0）t，迭代至第14次后的结果为当w二1时,x（14）二（3.0026645,3.9977796,-5.0005551”当w二1.24时，x（14）二（3.0000002,3.9999999,-5.）t可见取w二1.24时的收敛速度比w二1时的收敛速度要快。若要精确到小数后7位，当w二1（GS法）时需迭代31次，而当w二1

14、.24（SOR法）时只需迭代14次，它表明松弛因子w选取的好坏，对收敛速度影响很大。六程序设计及实验结果1.Jacobi迭代法，用mathematica编写程序如下：（i）计算迭代矩阵ClearEvaluateContext<>"*"A=5,0,0,-3,-1,-1,4,0,0,-1,0,0,2,-1,0,-1,0,0,4,-2,0,0,0,-1,2;b=2,3,-1,0,-1;L=TableIfi>j,Ai,j,0,i,LengthA,j,LengthTransposeA;D1=TableIfi=j,Ai,j,0,i,LengthA,j,LengthTr

15、ansposeA;U=TableIfi<j,Ai,j,0,i,LengthA,j,LengthTransposeA;I1=IdentityMatrixLengthA;BJ=I1-InverseD1.A/NfJ=InverseD1.b/N;运行结果为厂00.25000.60.2、0.25000B=0000.50J0.250000.5j0000.50丿(ii)计算方程组的解精确到小数点后7位时，迭代次数、最后一次迭代的结果、最后两次的相对误差。s=Table0,i,1,Lengthb,Table1,i,1,Lengthb;k=2;x1=x0=Table0,i,1,Lengthb;WhileN

16、ormsk-sk-1,2>10-6,x1=B.x0+f;x0=x1;s=AppendJJTos,x1;k+s=Deletes,1,2;n=LengthssnNormsn-sn-1,2/Normsn,2运行结果为迭代次数为n=32最后一次迭代结果为x(32)=(0.086957356,0.60869617,-0.65217336,-0.3043471,-0.65217336)t最后两次的相对误差IIX(32)一X(31)|/vcvc二6.4345694x10-7<10一IIx32II(iii)计算平均收敛速度以及渐进收敛速度u=MaxEigenvaluesBJ;R=-Logu/NRk

17、=Table-LogPowerNormMatrixPowerBJ,k,1/k/N,k,8；Rk=Table-LogPowerNormMatrixPowerBJ,k,1/k/N,k,8；date1=Tablei,Rki,LogRki,i,LengthRk;TableFormdatelATableHeadingsH1HAH2HAH3HAH4HAH5HAH6HAH7HiignjnRkHrHlnRkHrTableSpacing2,4rTableAlignments运行结果为渐进收敛速度为R=0.41301462平均收敛速度如表所示kRkInRK10.22314355-1.4999400020.2989

18、1850-1.2075843030.35473695-1.0363787040.35677909-1.0306385050.37710417-0.9752338260.37622879-0.9775578370.38679290-0.9498658680.38586107-0.95227789(iv)对平均收敛速度取对?牧后作最小一乘拟合图像date2=TableLogRki,LogRki+1,i,LengthRk-1;line=Fitdate2,1,t,tP=Dline,tp1=ListPlotdate2,PlotStylePointSize0.018;p2=Plotline,t,MinLo

19、gRk-10-1,MaxLogRk+10-1,PlotRange->MinLogRk-10-2,MaxLogRk+10-2;Showp1,p2运行结果如图所示2.GaussSeidel迭代法，用mathematica编写程序如下：(i)计算迭代矩阵ClearEvaluateContext<>"*"A=5,0,0,-3,-1,-1,4,0,0,-1,0,0,2,-1,0,-1,0,0,4,-2,0,0,0,-1,2;b=2,3,-1,0,-1;L=TableIfi>j,Ai,j,0,i,LengthA,j,LengthTransposeA;D1=Tab

20、leIfi=j,Ai,j,0,i,LengthA,j,LengthTransposeA;U=TableIfi<j,Ai,j,0,i,LengthA,j,LengthTransposeA;B=-InverseL+D1.U/NGf=InverseL+D1.b/N;G运行结果为0.2、0.300.550.275丿Poo0.60000.15B=0000.5G0000.15k0000.075(ii)计算方程组的解精确到小数点后7位时，迭代次数、最后一次迭代的结果、最后两次的相对误差。s=Table0,i,1,Lengthb,Table1,i,1,Lengthb;k=2;x1=x0=Table0,i

21、,1,Lengthb;WhileNormsk-sk-1,2>10-6,x1=B.x0+f;x0=x1;s=AppendGGTos,x1;k+s=Deletes,1,2;n=LengthssnNormsn-sn-1,2/Normsn,2运行结果为迭代次数为n=18最后一次迭代结果为x(18)=(0.086956842,0.60869579,-0.65217368,-0.30434763,-0.65217382)T最后两次的相对误差11X"X"11=5.5402987x10-7<10-6IIX18II(iii)计算平均收敛速度以及渐进收敛速度u=MaxEigenva

22、luesB;G-Logu/Ndate1=Tablei,Rki,LogRki,i,LengthRk;Rk=Table-;LogPowerNormMatrixPowerBGrkTableFormdatelATableHeadingsH1HAH2HAH3HAH4HAH5HAH6HAH7HiignjnRkHrHlnRkHrTableSpacing2,4rTableAlignments运行结果为渐进收敛速度为R=0.85566611平均收敛速度如表所示kRkInRK10.22314355-1.4999400020.35667494-1.0309304030.52300533-0.6481636240.6

23、0617053-0.5005939350.65606964-0.4214883360.68933572-0.3720268770.71309721-0.3381375380.73091832-0.31345356(iv)对平均收敛速度取对数后作最小二乘拟合图像如图所示date2=TableLogRki,LogRki+1,i,LengthRk-1line=Fitdate2,1,t,tP=Dline,tp1=ListPlotdate2,PlotStyleTPointSize0018;p2=Plotline,t,MinLogRk-10-1,MaxLogRk+10-1,PlotRange->Mi

24、nLogRk-10-2,MaxLogRk+10-2;Showp1,p23.SOR迭代法，用mathematica编写程序如下：(i)当w二1(GS法)时，计算其渐进收敛速度以及方程组的解精确到小数点后7位时的迭代次数、最后一次迭代的结果、最后两次的相对误差。ClearEvaluateContext<>"*"A=4,3,0,3,4,-1,0,-1,4;b=24,30,-24;L=TableIfi>j,Ai,j,0,i,LengthA,j,LengthTransposeA;D1=TableIfi=j,Ai,j,0,i,LengthA,j,LengthTrans

25、poseA;U=TableIfi<j,Ai,j,0,i,LengthA,j,LengthTransposeA;BG=-InverseL+D1.U/N;fG=InverseL+D1.b;RR=MaxEigenvaluesBG;-LogRR/NS=Table0,i,1,Lengthb,Table1,i,1,Lengthb;k=2;x1=x0=Table0,i,1,Lengthb;WhileNormsk-sk-1,2>10-6,x1=BG.x0+fG;x0=x1;s=AppendTos,x1;k+s=Deletes,1,2;n=LengthssnNormsn-sn-1,2/Normsn,

26、2运行结果为渐进收敛速度为R=0.47000363迭代次数为n=31最后一次迭代结果为x(31)=(3.0000009,3.9999992,-5.0000002)T最后两次的相对误差11x)x)11二1.0098429x10-7<10-6|x(31)|(ii)当w二1.24(SOR法)时，计算其渐进收敛速度以及方程组的解精确到小数点后7位时的迭代次数、最后一次迭代的结果、最后两次的相对误差。ClearEvaluateContext<>"*"A=4,3,0,3,4,-1,0,-1,4;b=24,30,-24;L=TableIfi>j,Ai,j,0,i,LengthA,j,LengthTransposeA;D1=TableIfi=j,Ai,j,0,i,LengthA,j,LengthTransposeA;U