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

研究解线性方程组迭代收敛速度一．实验目的科学研究与生产实践中许多问题都可归结为线性方程组的求解，高效求解线性方程组成为了许多科学与工程计算的核心．迭代法就是用某种极限过程去逼近线性方程组精确解的方法，该方法具有对计算机的存贮单元需求少，程序计算简单，原始系数矩阵在计算过程中不变等优点，是求解大型稀疏矩阵方程组的重要方法。常用的迭代法有Jacobi迭代法、Gauss—seidel迭代法、逐次超松驰法(SOR法)等。二．实验摘要由迭代法平均收敛速度与渐进收敛速度的关系引入近似估计法，即通过对迭代平均收敛速度取对数，然后利用Mathematica软件对其进行拟合，给出拟合函数，最终得到了Jacobi迭代法、Gauss—seidel法的平均收敛速度收敛到渐进收敛速度的近似收敛阶，以及逐次超松驰法(SOR法)的渐进收敛速度，且该法适用于其他迭代法收敛速度的估计。三．迭代法原理1．Jacobi迭代法(J法)设方程组bAx，其中，nnnnijRaA)(，。nRbx,A为可逆矩阵，可分裂为,UDLA其中，00001,21323121nnnnaaaaaaL0000,122311312nnnnaaaaaaUnnaaaD2211从而由bAx得到，bDxADIbDxADDbDxULDx111111)()()(令ADIBJ1,bDfJ1,由此可构造出迭代公式：JkJkfxBx)()1(令初始向量)0,...,0,0()0(x，即可得到迭代序列，从而逼近方程组的解这种方法称为Jacobi迭代法，其中JB称为Jacobi迭代矩阵。2.Gauss-Seidel迭代法(GS法)与Jacobi迭代法类似，将方程组bAx中的系数矩阵A分裂为,UDLA,其中ULD,,与前面相同。与Jacobi迭代法所不同的是，Gauss-Seidel迭代法将Jacobi迭代公式中的bUxLxDxkkk)()()1(改为bUxLxDxkkk)()1()1(从而bAx可写成矩阵形式bUxxDLkk)()1()(，若设1)(DL存在，则bDLUxDLxkk1)(1)1()()(，其中，UDLBG1)(,bDLf1)(，于是Gauss—Seidel迭代公式的矩阵形式为fxBxkGk)()1(。其中，GB称为式（1）的Gauss—Seidel迭代法的迭代矩阵。注：由于Gauss-Seidel迭代充分利用了迭代过程的新信息，一般来说，它的迭代效果要比Jacobi迭代好。3.逐次超松弛方法（SOR法）根据Gauss-Seidel迭代法的迭代原理)()()1(1)1(bUxLxDxkkk，我们将其修改为)1()()1()1(kkkxwxwx，即对)(kx和)1(kx加入一个权因子，在这里我们称w为迭代公式的松弛系数。令)()()1(1)1(bUxLxDxkkk（Gauss—Seidel迭代法），则)()1()1()()1(1)()1()()1(bUxLxwDxwxwxwxkkkkkk从而bwLDwxwUDwwLDxkk1)(1)1()(])1[()(])1[()(1wUDwwLDGw，1)(wLDwfw所以wkwkfxGx)()1(将其写成分量的形式我们称该公式为基于Gauss—Seidel迭代下的超松弛迭代公式。注：1)关于松弛系数w的选取，我们有以下性质：（i）.设方程组bAx，其中，nnnnijRaA)(，nRbx,，则解方程的SOR迭代法收敛的必要条件是20w；（ii）若A为正定的对称矩阵，且20w，则方程组bAx关于SOR迭代法是收敛的；（iii）若A为正定的对称三对角矩阵，JB为Jacobi迭代矩阵，若JB的谱半径1)(JB，记)(JB，则SOR法的最优松弛因子为2112bw且2104/)1(4)(2222)由松弛系数的性质可知，通常只有当方程组的系数矩阵A为正定的对称矩阵时我们才使用SOR法；3)当松弛系数1w时，SOR法记为GS方法；4)关于（iii）提到的谱半径定义为：假设A为n阶可逆矩阵，n,...,,21是它的n个特征值，我们称||max)(1iniA为A的谱半径。容易证明A的谱半径是A的所有诱导范数的下确界，即||||inf)(||||AA其中，矩阵的范数如：niijnjaA111||max||||，njijniaA11||max||||，)(||||2AAAT。下面我们就来讨论以上方法的迭代收敛速度，首先我们介绍收敛速度快慢的比较方法。四．迭代法收敛速度的比较1.这两种迭代方法收敛性与)(0kBk是否成立有关，且收敛速度与0kB的速度有关。当1)(B时，kB趋于零矩阵的速度有赖于)(B的大小。一般说来，)(B越小，则kB趋于零矩阵的速度越快，反之就越慢。通常，当1)(B时，可以用正数)(lnB的大小作为迭代法渐进收敛速度的度量。这时)(lnB越大，迭代法的收敛速度愈大。2.平均收敛速度与渐进收敛速度之间的联系对于收敛的迭代法...)2,1,0(1kfBxxkk，)||ln(||/1kkkBR称为平均收敛速度(它与所用的范数以及k有关)；)(lnBR称为渐进收敛速度。可以证明kkRRlim，因此当k时假设成立下列渐近关系式cRRpkk1(常数)，为估计平均收敛速度收敛到渐进收敛速度的阶，当k充分大时有如下近似形式：1,1kcRRkk两边取对数得：cRpRkklnlnln此式说明当k比较大时，1lnkR与kRln有近似的线性关系，而它们的斜率刚好为收敛阶P，这样可以通过当k比较大时作出1lnkR与kRln的拟合曲线来估计出P值。五．实验举例例1：考虑五阶方程组1202412342355454143521541xxxxxxxxxxxxx1.其Jacobi迭代法的迭代矩阵为05.00005.000025.005.000025.000025.02.06.0000JB渐进收敛速度为：41301462.0)66165261.0ln())(ln(JBR则迭代平均收敛速度kR（见表1）以及取对数后作最小二乘拟合图像（见图1）如下所示：表1kkRKRln10.22314355-1.4999400020.29891850-1.2075843030.35473695-1.0363787040.35677909-1.0306385050.37710417-0.9752338260.37622879-0.9775578370.38679290-0.9498658680.38586107-0.95227789图11.51.41.31.21.11.01.201.151.101.051.000.95即53345499.0)ln(44226443.0)ln(1kkRR从而得到收敛阶为p=0.442264432.该方程组的Gauss—Seidel迭代矩阵为：275.0075.000055.015.000005.00003.015.00002.06.0000GB其渐进收敛速度为85566611.0)425.0ln())(ln(GBR则迭代平均收敛速度kR(见表2)以及取对数后作最小二乘拟合图像(见图2)如下所示：表2kkRKRln10.22314355-1.4999400020.35667494-1.0309304030.52300533-0.6481636240.60617053-0.5005939350.65606964-0.4214883360.68933572-0.3720268770.71309721-0.3381375380.73091832-0.31345356图21.41.21.00.80.60.40.90.80.70.60.50.40.3即11593362.0)ln(58472143.0)ln(1kkRR从而得到收敛阶为p=0.58472143。注：在本例中，由于方程组的系数矩阵严格对角占优，故前述两种迭代过程均收敛，依实际迭代过程：对Jacobi迭代有67-32)31()32(10106.4345694||||||||xxx而Gauss—Seidel迭代有67-18)17()18(10105.5402987||||||||xxx，即二者均收敛，但后者更快一些。这与收敛阶p=0.5850.442之间关系一致。例2：考虑方程组244304324343232121xxxxxxx用SOR迭代公式可得取初试量为Tx)0,0,0()0(，迭代至第14次后的结果为当1w时，Tx)5.0005551-3.9977796,3.0026645,()14(当24.1w时，Tx)5.-3.9999999,3.0000002,()14(可见取24.1w时的收敛速度比1w时的收敛速度要快。若要精确到小数后7位，当1w(GS法)时需迭代31次，而当24.1w（SOR法）时只需迭代14次，它表明松弛因子w选取的好坏，对收敛速度影响很大。六．程序设计及实验结果1.Jacobi迭代法，用mathematica编写程序如下：（i）计算迭代矩阵Clear[Evaluate[Context[]*]]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=Table[If[ij,A[[i,j]],0],{i,Length[A]},{j,Length[Transpose[A]]}];D1=Table[If[i==j,A[[i,j]],0],{i,Length[A]},{j,Length[Transpose[A]]}];U=Table[If[ij,A[[i,j]],0],{i,Length[A]},{j,Length[Transpose[A]]}];I1=IdentityMatrix[Length[A]];BJ=I1-Inverse[D1].A//NfJ=Inverse[D1].b//N;运行结果为05.00005.000025.005.000025.000025.02.06.0000JB（ii）计算方程组的解精确到小数点后7位时，迭代次数、最后一次迭代的结果、最后两次的相对误差。s={Table[0,{i,1,Length[b]}],Table[1,{i,1,Length[b]}]};k=2;x1=x0=Table[0,{i,1,Length[b]}];While[Norm[s[[k]]-s[[k-1]],2]10-6,x1=BJ.x0+fJ;x0=x1;s=AppendTo[s,x1];k++]s=Delete[s,{{1},{2}}];n=Length[s]s[[n]]Norm[s[[n]]-s[[n-1]],2]/Norm[s[[n]],2]运行结果为迭代次数为n=32最后一次迭代结果为Tx)521733643471,-0.67336,-0.3017,-0.65216,0.6086960.08695735()32(最后两次的相对误差67-32)31()32(10106.4345694||||||||