Hilbert矩阵病态线性代数方程组的求解

实验一病态线性代数方程组的求解1.估计Hilbert矩阵2-条件数与阶数的关系运行tiaojianshu.m输入m=10可以得到如下表的结果2.选择不同维数，分别用Guass消去（LU分解），Jacobi迭代，GS迭代，SOR迭代求解，比较结果。说明：Hx=b，H矩阵可以由matlab直接给出，为了设定参考解，我们先设x为分量全1的向量，求出b,然后将H和b作为已知量，求x，与设定的参考解对比。对于Jacobi迭代，GS迭代，SOR迭代，取迭代初值x0为0向量，迭代精度eps=1.0e-6，迭代次数100000,SOR迭代中w=1.2和0.8分别计算。a.n=5x分量Gauss法J迭代GS迭代SOR迭代实际解xw=1.2w=0.8x(1)1.0000-8.81790.99980.99990.99981x(2)1.0000-Inf1.00281.00191.00311x(3)1.0000-Inf0.98790.99190.98631x(4)1.0000-Inf1.01811.01201.02091x(5)1.0000-Inf0.99120.99420.98971迭代次数14229221603147b.n=8x分量Gauss法J迭代GS迭代SOR迭代实际解xw=1.2w=0.8x(1)1.00001.18981.00011.00021.00001x(2)1.0000Inf0.99740.99610.99861x(3)1.0000Inf1.01361.02011.00821x(4)1.0000Inf0.97940.96820.98751x(5)1.0000Inf0.99821.00340.99691阶数12345条件数119.28524.051.55e+44.76e+5阶数678910条件数1.49e+74.75e+81.52e+104.93e+111.60e+13x(6)1.0000Inf1.01411.01631.01031x(7)1.0000Inf1.01011.01041.00891x(8)1.0000Inf0.98700.98510.98941迭代次数834278409473c.n=10x分量Gauss法J迭代GS迭代SOR迭代实际解xw=1.2w=0.8x(1)1.0000-0.50271.00011.00011.00021x(2)1.0000-1.17630.99820.99880.99741x(3)1.0000-1.64651.00561.00291.00871x(4)1.0000-Inf0.99931.00310.99611x(5)0.9999-Inf0.99160.99080.99041x(6)1.0004-Inf0.99680.99640.99741x(7)0.9994-Inf1.00521.00441.00681x(8)1.0006-Inf1.00871.00811.01021x(9)0.9997-Inf1.00391.00391.00411x(10)1.0001-Inf0.99040.99140.98861迭代次数269512960821276d.n=15x分量Gauss法J迭代GS迭代SOR迭代实际解xw=1.2w=0.8x(1)1.00001.60600.99990.99980.99991x(2)1.0000Inf1.00351.00381.00231x(3)0.9995Inf0.98210.98190.98771x(4)1.0059Inf1.02361.02361.01531x(5)0.9656Inf1.00741.00431.00821x(6)1.0813Inf0.99020.99270.99361x(7)1.1155Inf0.98560.98850.98751x(8)-0.3001Inf0.99050.99190.99051x(9)4.8294Inf0.99900.99850.99811x(10)-4.9705Inf1.00681.00491.00561x(11)6.0398Inf1.01131.00891.01041x(12)-0.5735Inf1.01131.00931.01081x(13)0.1389Inf1.00661.00581.00661x(14)1.8886Inf0.99750.99840.99791x(15)0.7794Inf0.98440.98740.98521迭代次数144982877515243取不同的n值，得到如下结果:对于Guass法，可以看出来，随着n的增大，求解结果误差变大，这是因为随着n增大，系数矩阵的条件数变大，微小的扰动就容易造成很大的误差。最后得不到精确解。对于Jacobi迭代，计算结果为Inf，说明是发散的。对于GS迭代和SOR迭代，结果是收敛的，但是可以看出迭代次数比较多，并且对于不同维数GS和SOR收敛速度不一样，有时候GS快，有时SOR快。对SOR取不同的w迭代速度也不一样，存在一个最优的松弛因子w。并且可以知道，迭代次数多少跟初值x0也有关系。3.讨论病态问题求解的算法。通过上面的实验分析，可以看出，求解病态矩阵的时候要小心，否则可能得不到所要求的精确度。可以采用高精度运算，用双倍多倍字长，使得由于误差放大而损失若干有效数字位之后，还能保留一些有效位。另外可以通过对原方程作某些预处理，降低系数矩阵的条件数，因为cond(aA)=cond(A),所以不能通过将每一个方程乘上相同的常数来达到这个目标，可考虑将矩阵的每一行和每一列分别乘上不同的常数，亦即找到可逆的对角阵D1和D2将方程组化为D1AD2y=D1b，x=D2y这称为矩阵的平衡问题，但是这样计算量比原问题本身要多。或者通过变分原理将求解线性方程组的问题转化为等价的求解无约束函数最优化问题的极小值等等，可以参考[1]郑洲顺，黄光辉，杨晓辉求解病态线性方程组的混合算法实验一所编程序如下：1.求条件数tiaojianshu.mm=input('inputm:=');N=[1:m];fori=1:mn=N(i);H=hilb(n);k=cond(H);disp('矩阵的阶数')disp(n)disp('矩阵')disp(H)disp('矩阵的条件数')disp(k)end2.Guass法①Guass.mfunctionguass(n)n=input('请输入系数矩阵的维数n=');H=hilb(n);x_exact=ones(n,1);b=H*x_exact;x=Doolittle(H,b)a=input('是否继续,继续请按1,结束请按0:');ifa==1guass(n);elseend②Doolittle.mfunctionx=Doolittle(A,b)%LU分解求Ax=b的解N=size(A);n=N(1);L=eye(n,n);%L的对角元素为1U=zeros(n,n);%U为零矩阵U(1,1:n)=A(1,1:n)%U第一行L(1:n,1)=A(1:n,1)/U(1,1)%L第一列fork=2:nfori=k:nU(k,i)=A(k,i)-L(k,1:(k-1))*U(1:(k-1),i)endforj=(k+1):nL(j,k)=(A(j,k)-L(j,1:(k-1))*U(1:(k-1),k))/U(k,k)endendy=solveDownTriangle(L,b);%调用下三角系数矩阵求解函数x=solveUpTriangle(U,y);%调用上三角系数矩阵求解函数③solveDownTriangle.mfunctionx=solveDownTriangle(A,b)%求下三角系数矩阵的线性方程组Ax=bN=size(A);n=N(1);fori=1:nif(i1)s=A(i,1:(i-1))*x(1:(i-1),1);elses=0;endx(i,1)=(b(i)-s)/A(i,i);end④solveUpTriangle.mfunctionx=solveUpTriangle(A,b)%求上三角系数矩阵的线性方程组Ax=bN=size(A);n=N(1);fori=n:-1:1if(in)s=A(i,(i+1):n)*x((i+1):n,1);elses=0;endx(i,1)=(b(i)-s)/A(i,i);end3.Jacobi法function[x,n]=jacobi(a,x0)a=input('请输入系数矩阵的维数n=');x0=input('请输入迭代初始向量x0');eps=1.0e-6;%解的精度控制m=10000;%迭代步数控制H=hilb(a);%生成h矩阵x_exact=ones(a,1);%求出x精确值b=H*x_exact;D=diag(diag(H));L=-tril(H,-1);U=-triu(H,1);B=D\(L+U);f=D\b;x=x0;n=0;tol=1;whiletol=epsx=B*x0+f;n=n+1;tol=norm(x-x0);x0=x;if(n=m)disp('迭代次数太多，可能不收敛');return;endenddisp(x)disp(n)4.GS法function[x,n]=gauseidel(a,x0)a=input('请输入系数矩阵的维数n=');x0=input('请输入迭代初始向量x0');eps=1.0e-6;%解的精度控制m=100000;%迭代步数控制H=hilb(a);%生成h矩阵x_exact=ones(a,1);%求出x精确值b=H*x_exact;D=diag(diag(H));L=-tril(H,-1);U=-triu(H,1);G=(D-L)\U;f=(D-L)\b;x=x0;n=0;tol=1;whiletol=epsx=G*x0+f;n=n+1;tol=norm(x-x0);x0=x;if(n=m)disp('迭代次数太多，可能不收敛');return;endenddisp(x)disp(n)5.SOR法function[x,n]=SOR(a,x0,w)a=input('请输入系数矩阵的维数n=');x0=input('请输入迭代初始向量x0=');w=input('请输入w=');eps=1.0e-6;%解的精度控制m=100000;%迭代步数控制H=hilb(a);%生成h矩阵x_exact=ones(a,1);%求出x精确值b=H*x_exact;if(w=0||w=2)error;return;endD=diag(diag(H));L=-tril(H,-1);U=-triu(H,1);B=inv(D-L*w)*((1-w)*D+w*U);f=w*inv((D-L*w))*b;x=x0;n=0;tol=1;whiletol=epsx=B*x0+f;n=n+1;tol=norm(x-x0);x0=x;if(n=m)disp('迭代次数太多，可能不收敛');return;endenddisp(x)disp(n)