病态线性方程的求解

（Hilbert 矩阵）病态线性方程组的求解理论分析表明，数值求解病态线性方程组很困难。

考虑求解如下的线性方程组的求解Hx = b ，期中H 是Hilbert 矩阵，()ij n n Hh ，11ij h i j ，i ，j = 1,2，…，n 1.估计矩阵的2条件数和阶数的关系2.对不同的n ，取(1,1,,1)nx K ?，分别用Gauss 消去，Jacobi 迭代，Gauss-seidel 迭代，SOR 迭代和共轭梯度法求解，比较结果。

3.结合计算结果，试讨论病态线性方程组的求解。

第1小题：condition.m %第1小题程序t1=20;%阶数n=20x1=1:t1;y1=1:t1;for i=1:t1H=hilb(i);y1(i)=log(cond(H));endplot(x1,y1);xlabel('阶数n');ylabel('2-条件数的对数（log(cond(H)）');title('2-条件数的对数（log(cond(H)）与阶数n 的关系图');t2=200;%阶数n=200x2=1:t2;y2=1:t2;for i=1:t2H=hilb(i);y2(i)=log(cond(H));endplot(x2,y2);xlabel('阶数n');ylabel('2-条件数的对数（log(cond(H)）');title('2-条件数的对数（log(cond(H)）与阶数n 的关系图');画出Hilbert 矩阵2-条件数的对数和阶数的关系n=200时n=20时从图中可以看出，1）在n小于等于13之前，图像近似直线log（cond（H））~1.519n-1.8332）在n大于13之后，图像趋于平缓，并在一定范围内上下波动，同时随着n的增加稍有上升的趋势第2小题：solve.m%m第2小题主程序N=4000;xGauss=zeros(N,1);xJacobi=zeros(N,1);xnJ=zeros(N,1);xGS=zeros(N,1);xnGS=zeros(N,1);xSOR=zeros(N,1);xnSOR=zeros(N,1);xCG=zeros(N,1);xnCG=zeros(N,1);for n=1:N;x=ones(n,1);t=1.1;%初始值偏差x0=t*x;%迭代初始值e=1.0e-8;%给定的误差A=hilb(n);b=A*x;max=100000000000;%可能最大的迭代次数w=0.5;%SOR迭代的松弛因子G=Gauss(A,b);[J,nJ]=Jacobi(A,b,x0,e,max);[GS,nGS]=G_S(A,b,x0,e,max);[S_R,nS_R]=SOR(A,b,x0,e,max,w);[C_G,nC_G]=CG(A,b,x0,e,max);normG=norm(G'-x);xGauss(n)=normG;normJ=norm(J-x);nJ;xJacobi(n)=normJ;xnJ(n)=nJ;normGS=norm(GS-x);nGS;xGS(n)=normGS;xnGS(n)=nGS;normS_R=norm(S_R-x);nS_R;xSOR(n)=normS_R;xnSOR(n)=nS_R;normC_G=norm(C_G-x);nC_G;xCG(n)=normC_G;xnCG(n)=nC_G;endGauss.m%Gauss消去法function x=Gauss(A,b)n=length(b);l=zeros(n,n);x=zeros(1,n);%消去过程for i=1:n-1for j=i+1:nl(j,i)=A(j,i)/A(i,i);for k=i:nA(j,k)=A(j,k)-l(j,i)*A(i,k);endb(j)=b(j)-l(j,i)*b(i);endend%回代过程x(n)=b(n)/A(n,n);for i=n-1:-1:1c=A(i,:).*x;x(i)=(b(i)-sum(c(i+1:n)))/A(i,i);endJacobi.m%Jacobi迭代,x0表示迭代初值,e表示允许误差(迭代停止条件),n表示迭代次数,m 可能最大的迭代次数function [x,n]=Jacobi(A,b,x0,e,m)n=length(A);D=diag(diag(A));U=-triu(A,1);L=-tril(A,-1);B=D\(L+U);f=D\b;x=B*x0+f;n=1;while norm(x-x0)>ex0=x;x=B*x0+f;n=n+1;if n>mdisp('Jacobi迭代次数过多，迭代可能不收敛');break;endendG_S.m%Gauss-Seidel迭代,x0表示迭代初值,e表示允许误差(迭代停止条件),n表示迭代次数,m可能最大的迭代次数function [x,n]=G_S(A,b,x0,e,m)n=length(A);D=diag(diag(A));U=-triu(A,1);L=-tril(A,-1);B=(D-L)\U;f=(D-L)\b;x=B*x0+f;n=1;while norm(x-x0)>ex0=x;x=B*x0+f;n=n+1;if n>mdisp('Gauss-Seidel迭代次数过多，迭代可能不收敛');break;endendSOR.m%SOR超松弛迭代,x0表示迭代初值,e表示允许误差(迭代停止条件),n表示迭代次数,m可能最大的迭代次数,w松弛因子function [x,n]=SOR(A,b,x0,e,m,w)n=length(A);D=diag(diag(A));U=-triu(A,1);L=-tril(A,-1);B=(D-w*L)\((1-w)*D+w*U);f=(D-w*L)\b*w;x=B*x0+f;n=1;while norm(x-x0)>ex0=x;x=B*x0+f;n=n+1;if n>mdisp('SOR超松弛迭代次数过多，迭代可能不收敛');break;endendCG.m%CG共轭梯度法,x0表示迭代初值,e表示允许误差(迭代停止条件),n表示迭代次数,m可能最大的迭代次数function [x,n]=CG(A,b,x0,e,m)r=b-A*x0;p=r;alpha=(r'*r)/(p'*(A*p));x=x0+alpha*p;r1=b-A*x;n=1;while norm(r1)>ebelta=(r1'*r1)/(r'*r);p=r1+belta*p;r=r1;x0=x;alpha=(r'*r)/(p'*(A*p));x=x0+alpha*p;r1=b-A*x;n=n+1;if n>mdisp('CG共轭梯度法迭代次数过多，迭代可能不收敛');break;endend。

第2次作业病态线性方程组的计算在科学技术、工程和经济等领域中都会遇到求解线性方程组问题。

在很多有广泛背景的数学问题中，如样条逼近、微分方程边值问题的差分法和有限元方法都要求求解线性方程组，而且通常会遇到求解大型方程组的问题。

在数值求解大型方程的过程中，一般会遇到很到计算上的问题，矩阵的病态性是其中之一。

矩阵的病态性质一般来自两个方面，一是由于矩阵本身的性质所引起，如希尔伯特矩阵。

这种矩阵，虽然阶数可能不大，但其条件数却非常大；另一类病态性质是由于矩阵的阶数所引起的。

矩阵在阶数比较小时，性质很好；但当阶数升高后，条件数则变得很大。

如用差分法或有限元法离散微分方程后得到的线性方程组。

它们的理论性质很好，甚至在理论上可证明它们是对称正定的。

然而，由于矩阵得阶数很高，可达到上百万阶，因而矩阵得条件数是相当大的。

方程的病态性质给方程组的数值求解带来巨大麻烦，甚至使某些看似理论性质很好的方程如对称正定矩阵，也不太可能得到通常意义下唯一的数值解。

病态矩阵是指矩阵有很小的扰动时，所引起解的变化很大的矩阵。

矩阵的病态性质一般可用矩阵的条件数来指示。

矩阵的条件数记为：-(cond)=AAA虽然矩阵的病态性，如上面提到的，可能来自两个方面,但在与具体问题相联系的求解过程中遇到的问题是基本一致。

下面这个方程组来自对偏微分方程（Stokes方程）的某种有限元方法离散后得到的。

该微分方程为⎪⎩⎪⎨⎧=Ω∈=Ω∈=∇+∆Ω∂0),( , 0 div ),( , ),(u y x u y x y x f p u 其中)),(),,((),(21y x u y x u y x u =为二维空间的向量函数，为标量函数),(y x p 。

为方便起见，现已经有方程组的系数矩阵和右端向量：这些数据存于数据文件‘stokes_coff_....mat ’中，可用Matlab 的load 命令调入即可。

关于数据文件中数据的说明：CoffA:方程组的系数矩阵；具有如下形式：⎪⎪⎭⎫ ⎝⎛=C B B ACoffA T变量BT 即为T B ；CoffA 为对称正定矩阵。

病态线性代数方程组的求解理论的分析表明，求解病态的线性代数方程组是困难的。

考虑方程组Hx = b 的求解，其中H 为Hilbert 矩阵，n n ij h H ⨯=)(，11-+=j i h ij ，n j i ,...,2,1,=1. 估计Hilbert 矩阵2-条件数与阶数的关系；2. 选择问题的不同维数，分别用Gauss 消去法，Jacobi 迭代，GS 迭代和SOR 迭代求解，比较结果；3. 讨论病态问题求解的算法。

解： 1、取Hilbert 矩阵阶数最高分别为n=20和n=100。

采用Hilbert 矩阵的2-条件数作为实验的比较对象，画出的曲线如下图所示：lg(())n cond H n从图中可以看出，在n ≤13之前，图像接近直线，在n ＞13之后，图像趋于平缓，在一定的范围内上下波动。

为了比较图像的线性部分，作出一条直线与已知曲线进行比较。

比较直线的关系式为：833.1519.1))(lg(-=n H cond n ，结果下图所示。

nl g (c o n d (H n ))lg(cond(Hn))~n 关系图从图2中可以看出，当n 较小时，n H cond n ~))(lg(之间近似满足线性关系。

当n 继续增大到100时，n H cond n ~))(lg(关系图下图所示：从图中可以看出，图像的走势符合在n=20时的猜想，在n 大于一定的值之后，图像趋于平缓，且在一定范围内震荡，同时又有一定上升趋势，但上升速度很慢。

2、选择不同的阶数n ，设方程组的精确解为xz=(1,1,…,1)T进行计算，用四种方法解x_Guass1、x_Jacobi1、x_GS1、x_SOR1对比表如下nl g (c o n d (H n ))lg(cond(Hn))~n关系图nl g (c o n d (H n ))lg(cond(Hn))~n 关系图Gauss消去法求解：选择问题的阶数为3~8时，用Gauss消去法求得的解与精确解一致，当阶数为9~14时，解开始出现偏差，而且n越大，偏差越大。

实验一病态线性代数方程组的求解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分别计算。

取不同的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 ('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('请输入系数矩阵的维数n=');H=hilb(n);x_exact=ones(n,1);b=H*x_exact;x=Doolittle(H,b)a=input('是否继续,继续请按1,结束请按0:');if a==1guass(n);else end②Doolittle.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:nfor i=k:nU(k,i)=A(k,i)-L(k,1:(k-1))*U(1:(k-1),i)endfor j=(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.mfunction x=solveDownTriangle(A,b)%求下三角系数矩阵的线性方程组Ax=bN=size(A);n=N(1);for i=1:nif (i>1)s=A(i,1:(i-1))*x(1:(i-1),1);elses=0;endx(i,1)=(b(i)-s)/A(i,i);end④solveUpTriangle.mfunction x=solveUpTriangle(A,b)% 求上三角系数矩阵的线性方程组Ax=bN=size(A);n=N(1);for i=n:-1:1if (i<n)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;while tol>=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;while tol>=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;while tol>=epsx=B*x0+f;n=n+1;tol=norm(x-x0);x0=x;if (n>=m)disp('迭代次数太多，可能不收敛');return;endenddisp(x)disp(n)。

数值实验3_3 病态的线性方程组求解
自动化系李琳琳 2004211068
1．选择问题的维数为6，分别用Gauss消去法、J迭代方法、GS迭代方法和SOR 迭代方法求解并与问题的解比较，所得结果见下表
2．逐步增大问题的维数，仍用上述方法求解，结果如下：
由上述计算结果可以看出：
病态方程组的数值求解必须小心进行，否则得不到所要求的准确度或不稳定。

由本题中所做的数值实验可以看出，对于系数矩阵为Hilbert矩阵的病态方程组，Guass(即LU分解)方法和J迭代法都是无效的，结果发散。

而GS迭代和SOR迭代方法则较有效的解决了这个问题，得到较为精确的结果。

另一方面，也可以说明GS方法和SOR方法在收敛性方面更有优势。

其中，当系数矩阵A对称正定时，GS法一定收敛，而J法却不一定；且采用最优松弛因子的SOR方法要比GS法和J法收敛快得多。

安徽工业大学数理科学与工程学院病态线性方程组的求解专业数学与应用数学班级数***班学号*******姓名 ***指导教师***二○一五年五月一、设计目的：为了更加透彻的熟悉数值分析课程，学习各种数学软件的使用，锻炼自己对知识的实际运用能力。

二、引言：用直接法求解AX=F 线性方程组,对于系数矩阵A 对角占优是很有效的。

方程阶数不高时,人们经常使用;而当方程组阶数大时,由于积累误差,导致结果失真。

为了克服误差积累问题,通常用迭代法。

它具有可达到所要求的精度和对计算内存要求不大的优点,对求解大型线性方程组,迭代法计算时间远比直接法少,所以在实际计算中,迭代法也被人们广泛使用。

然而迭代法要研究迭代格式的收敛性,如Jacobi 迭代对系数矩阵为病态矩阵不收敛,为此我们提供一种修改的Jacobi 迭代,并给出一些数值例子来说明有较好的效果。

三、解线性方程组迭代法的描述设线性方程组AX=P ，这里A:{a ij },X:{x i },F:{f i },1≤i,j ≤n,为了更广泛地应用，对A 只限制实的非奇异矩阵，那么，若给定初值x )0(,我们熟知的有: Jacobi 迭代:x )1(+k i=(f i -∑njk j ij x a )()/a ii j i ≠, 1≤i ≤n四、求解病态线性方程组的另一种迭代解法设线性方程BX=F, 这里系数矩阵B 是病态的，指的是矩阵条件数是较大的。

条件数越大，就越难求得准确解，为此，我们将方程的两端同加DX 项，那 么相应的Jacobi 迭代有：X ])([)()(1)1(k k X H D F A D -++=-+ (1)这里，A 为B 的对角阵，即A: {b ii },H:{b ij } j i ≠ 1≤i,j ≤n,记M= )()(1H D A D -+-，)()1()(k k k X X X -=∆+,那么有：)0()1()(X M X M X k k k ∆==∆=∆-Λ由此可见（1）式收敛，迭代矩阵M 的谱半径应满足p(M )<1，谱半径若用M 的特征值判断，则较为繁锁；若B 不可约，根据线性代数对角占优简单迭代必收敛的性质，为简单起见，我们取D 为对角线阵，即D: {d i },为保证收敛就得取d i =Sign{∑jii ij b b |,|}，符号Sign{a,b}的含义是与b 同号，数值取a,这是充分条件，实际计算中有时可放宽处理，比如可取d i =Sign{ii ij ijb b Max |,|}或者d i =Sign{ii ij jb b Max |,|} j i ≠，因而相应的Jacobi 迭代修改为：x )1(+k i =(f i -∑≠+ij k i i k j ij x d x b )()()/(b ii +d i ) （2）下面列举普通Jacobi 迭代不收敛，解不出正确结果，而用修改的Jacobi 迭代可求出满意结果的例子：例 1：A=⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡--211111112 F=⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡--032 精确解为X=[]T 111---，普通Jacobi 迭代不收敛，取d i =Sign{ii ij ja a Max |,|} j i ≠， 1≤i,j ≤3用（2）式迭代40步，X=[]T000003.1000053.1000035.1---例 2：方程系数为Hilbert 阵H:{)1/(1-+=j i h ij } F:{∑=jij i h f } 1≤i,j ≤n精确解X=[]T 111Λ，普通Jacobi 迭代发散，取d i =Sign{∑jii ij h h |,|}用修改的（2）式迭代，n=1200时迭代10870步结果摘录如下：五、结论与问题以上数值试验表明该迭代算法对求解病态线性方程组是有效的，其优点是可达到预定的精度。

数值分析课程实验报告 实验名称 病态线性方程组的算法设计
班级
学号 姓名 序号
任课教师 评分 一、 实验目的
1、 初步病态线性方程组的判定。

2、 初步了解常规方法在求解病态线性方程组时遇到的困难。

3、 针对病态问题设计求解算法并验证算法的有效性。

二、用文字或图表记录实验过程和结果
1、Hilbert 矩阵如下：
11/21/1/21/31/(1)()1/1/(1)1/(21)n ij n
n H h n n n ⎡⎤⎢⎥+⎢⎥==⎢⎥⎢⎥+-⎣⎦L L M M M L 其中1(1)ij h i j =+-，它是一个对称正定矩阵，并且()n cond H 随着n 的增加迅速增加，利用Matlab 分析如下：
可以发现在阶数不断增大
Hilbert 矩阵的条件数不断增大，
这样使得求解Hilbert 病态方程
变得非常困难，即使A 或b 有微小
扰动时，即使求解过程是精确进
行的（即没有舍入误差），所得的
解相对于原方程的解也会有很大
的相对误差。

这就需要提出病态
线性方程组的求解方法，对于一
般的方程求解常用的有高斯（选
主元）消去算法、高斯—赛德尔迭代。

本试验先使用用列主元高斯消去算法和高斯-赛德尔迭代算法求解线性方程组：
n
H x b = 其中11(,,),(1,2,,)n T n i ij j b b b b h i n ====∑L L 。

2、高斯列主元消去算法
(1)设计流程图：。