《Hilbert矩阵病态线性代数方程组的求解.docx》由会员分享，可在线阅读，更多相关《Hilbert矩阵病态线性代数方程组的求解.docx（8页珍藏版）》请在三一文库上搜索。
1、.实验一 病态线性代数方程组的求解1.估计Hilbert矩阵2-条件数与阶数的关系运行tiaojianshu.m 输入m=10 可以得到如下表的结果阶数12 345条件数119.28524.051.55e+44.76e+5阶数678910条件数1.49e+74.75e+81.52e+104.93e+111.60e+132.选择不同维数，分别用Guass消去（LU分解），Jacobi迭代，GS迭代，SOR迭代求解，比较结果。说明：Hx=b，H矩阵可以由matlab直接给出，为了设定参考解，我们先设x为分量全1的向量，求出b,然后将H和b作为已知量，求x，与设定的参考解对比。对于Jacobi迭代，
2、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迭代次数142292216031
3、47b. 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.99691x(6)1.0000Inf1.01411.01631.01031x(7)1.0000Inf1.01011.01041.00891x(8)1.0000Inf0.98700.98510.989
4、41迭代次数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.00681
5、x(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.96
6、56Inf1.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.8886In
7、f0.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。并且可以知道，迭代
8、次数多少跟初值x0也有关系。3.讨论病态问题求解的算法。 通过上面的实验分析，可以看出，求解病态矩阵的时候要小心，否则可能得不到所要求的精确度。可以采用高精度运算，用双倍多倍字长，使得由于误差放大而损失若干有效数字位之后，还能保留一些有效位。另外可以通过对原方程作某些预处理，降低系数矩阵的条件数，因为cond(aA)=cond(A),所以不能通过将每一个方程乘上相同的常数来达到这个目标，可考虑将矩阵的每一行和每一列分别乘上不同的常数，亦即找到可逆的对角阵D1和D2将方程组化为D1AD2y=D1b，x=D2y这称为矩阵的平衡问题，但是这样计算量比原问题本身要多。或者通过变分原理将求解线性方程组的
9、问题转化为等价的求解无约束函数最优化问题的极小值等等，可以参考1郑洲顺，黄光辉，杨晓辉 求解病态线性方程组的混合算法实验一所编程序如下：1.求条件数tiaojianshu.mm=input ('input m:=') ;N=1:m;for i=1:mn=N(i);H=hilb(n);k=cond(H);disp('矩阵的阶数')disp(n)disp('矩阵')disp(H)disp('矩阵的条件数')disp(k)end2.Guass法Guass.mfunction guass(n)n=input('请输入系数矩阵的维数
10、n=');H=hilb(n);x_exact=ones(n,1);b=H*x_exact;x=Doolittle(H,b)a=input('是否继续,继续请按1,结束请按0:');if a=1 guass(n);else endDoolittle.mfunction x=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第一列for k=2:n for i=
11、k:n U(k,i)=A(k,i)-L(k,1:(k-1)*U(1:(k-1),i) end for j=(k+1):n L(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.mfunction x=solveDownTriangle(A,b)%求下三角系数矩阵的线性方程组Ax=bN=size(A);n=N(1);for i=1:n if (i>1)
12、s=A(i,1:(i-1)*x(1:(i-1),1); else s=0; end x(i,1)=(b(i)-s)/A(i,i);endsolveUpTriangle.mfunction x=solveUpTriangle(A,b)% 求上三角系数矩阵的线性方程组Ax=bN=size(A);n=N(1);for i=n:-1:1 if (i<n) s=A(i,(i+1):n)*x(i+1):n,1); else s=0; end x(i,1)=(b(i)-s)/A(i,i);end3.Jacobi法function x,n=jacobi(a,x0)a=input('请输入系数矩阵
13、的维数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=Db;x=x0;n=0;tol=1;while tol>=eps x=B*x0+f; n=n+1; tol=norm(x-x0); x0=x; if (n>=m) disp('迭代次数太多，可能不收敛');
14、 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;while tol>
15、;=eps x=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;while tol>=eps x=B*x0+f; n=n+1; tol=norm(x-x0); x0=x; if (n>=m) disp('迭代次数太多，可能不收敛'); return; endenddisp(x)disp(n)*;