数值实验 病态的线性方程组的求解

实验五 解线性方程组的直接方法实验5.1 （主元的选取与算法的稳定性）问题提出：Gauss 消去法是我们在线性代数中已经熟悉的。

但由于计算机的数值运算是在一个有限的浮点数集合上进行的，如何才能确保Gauss 消去法作为数值算法的稳定性呢？Gauss 消去法从理论算法到数值算法，其关键是主元的选择。

主元的选择从数学理论上看起来平凡，它却是数值分析中十分典型的问题。

实验内容：考虑线性方程组n n n R b R A b Ax ∈∈=⨯,,编制一个能自动选取主元，又能手动选取主元的求解线性方程组的Gauss消去过程。

实验要求：（1）取矩阵⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡=⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡=1415157,6816816816M O O O b A ，则方程有解T x )1,,1,1(*Λ=。

取n=10计算矩阵的条件数。

让程序自动选取主元，结果如何？（2）现选择程序中手动选取主元的功能。

每步消去过程总选取按模最小或按模尽可能小的元素作为主元，观察并记录计算结果。

若每步消去过程总选取按模最大的元素作为主元，结果又如何？分析实验的结果。

（3）取矩阵阶数n=20或者更大，重复上述实验过程，观察记录并分析不同的问题及消去过程中选择不同的主元时计算结果的差异，说明主元素的选取在消去过程中的作用。

（4）选取其他你感兴趣的问题或者随机生成矩阵，计算其条件数。

重复上述实验，观察记录并分析实验结果。

思考题一：（Vadermonde 矩阵）设⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡=⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡=∑∑∑∑====n i i n n i i n i i n i i n n n n n n n x x x x b x x x x x x x x x x x x A 002010022222121102001111M ΛM ΛM M M ΛΛΛ，， 其中，n k k x k ,,1,0,1.01Λ=+=，（1）对n=2,5,8，计算A 的条件数；随n 增大，矩阵性态如何变化？（2）对n=5，解方程组Ax=b ；设A 的最后一个元素有扰动10-4，再求解Ax=b（3）计算（2）扰动相对误差与解的相对偏差，分析它们与条件数的关系。

第一次作业——病态问题实验1.1病态问题实验目的：算法有“优”和“劣”之分，问题也有“好”和“坏”之别，对数值方法的研究而言，所谓坏问题就是问题本身对扰动敏感，反正属于好问题。

问题提出：考虑一个高次的代数多项式∏=-=-⋅⋅⋅--=201) ()20()2)(1()(kk xxxxxp显然该多项式的全部根为1，2，……，20，共计20个，且每个根都是单重的。

现考虑该多项式的一个扰动)(19=+xxpε，其中，ε是一个非常小的数，这相当于是对方程中的x19的系数做一个小的扰动，比较两个方程根的差别，从而分析方程的解对扰动的敏感性。

function t_charp1_1%数值实验1.1变态问题%输入：[0 20]之间的扰动项及小的扰动常数%输出：加扰动都得到的全部根clcresult=inputdlg({'请输入扰动项：在[0 20]之间的整数：'},'charpt 1-1',1,{'19'});Numb=str2num(char(result));if((Numb>20)|(Numb<0)) errordlg('请输入正确的扰动项:[0 20]之间的整数!');return;endresult=inputdlg({'请输入（0 1）间的扰动常数：'},'charpt 1_1',1,{'0.00001'});ess=str2num(char(result));ve=zeros(1,21);ve(21-Numb)=ess;root=roots(poly(1:20)+ve);disp(['对扰动项',num2str(Numb),'加扰动',num2str(ess),'得到的全部根为：']);disp(num2str(root));结果分析：请输入扰动项：在[0 20]之间的整数：19 请输入（0 1）间的扰动常数：0.00001 对扰动项19加扰动1e-005得到的全部根为：22.5961+2.3083i22.5961-2.3083i18.8972+5.00563i18.8972-5.00563i14.9123+4.95848i14.9123-4.95848i12.0289+3.73551i12.0289-3.73551i10.059+2.33019i10.059-2.33019i8.63833+1.05643i8.63833-1.05643i7.70881+0i7.02809+0i5.99941+0i5.00001+0i4+0i3+0i 2+0i1+0i请输入扰动项：在[0 20]之间的整数：14 请输入（0 1）间的扰动常数：0.00001 对扰动项14加扰动1e-005得到的全部根为：2019.000317.99817.006215.987915.014913.988513.004612.000910.997610.00178.999328.000176.99997654321请输入扰动项：在[0 20]之间的整数：13 请输入（0 1）间的扰动常数：0.00001 对扰动项13加扰动1e-005得到的全部根为：20.000218.998718.003816.995215.999515.011713.979113.020711.986811.00589.998469.000188.000046.999986543 21请输入扰动项：在[0 20]之间的整数：1 请输入（0 1）间的扰动常数：0.00001 对扰动项1加扰动1e-005得到的全部根为：20.000218.998718.003816.995215.999515.011713.979113.020711.986811.00589.998469.000188.000046.99998654321结果分析：如图1、2所示，由上述结果可以看出，对加相同的干扰值ε在X的不同项上与不加干扰得到的结果相比可以看出，对扰动项x19加扰动1e-005得到的全部根、对扰动项x14加扰动1e-005得到的全部根与没有加扰动时的结果相比，相差较大，而且对扰动项x19加1e-005的结果干扰较大。

（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矩阵，H=(hij)n⨯n，hij=1，i,j=1,2,...,n i+j-11. 估计Hilbert矩阵2-条件数与阶数的关系；2. 选择问题的不同维数，分别用Gauss消去法，Jacobi迭代，GS迭代和SOR迭代求解，比较结果；3. 讨论病态问题求解的算法。

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

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

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

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

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

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

2、选择不同的阶数n，设方程组的精确解为xz=(1,1,…,1)T进行计算，用四种方法解x_Guass1、x_Jacobi1、x_GS1、x_SOR1对比表如下表所示。

Gauss消去法求解：选择问题的阶数为3~8时，用Gauss消去法求得的解与精确解一致，当阶数为9~14时，解开始出现偏差，而且n越大，偏差越大。

用迭代法求解：取初始向量为1.2(1,1,…,1)T.无论n为多少阶，用Jacobi迭代方法迭代出现发散的不稳定现象，无法求解；用GS迭代方法迭代不发散，能求得解，但收敛非常缓慢，当迭代次数取得相当大（20000次）时解仍在精确解附近波动；取w=1.5，用SOR迭代方法迭代不发散，能求得解，收敛速度较GS迭代快一些，但仍非常缓慢。

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

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

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

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

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

数值分析实验课第二单元第二部分数据拟合的数值方法实验报告问题重述：发现并观察数据拟合法正规方程组的病态问题.用拟合法求a ,b 的方法如下：⎪⎪⎩⎪⎪⎨⎧=+=+∑∑∑∑∑=====ni i i n i i n i i ni i n i i y x b x a x y b x na 112111其中n 为数据的个数 例如课本p50用拟合法求a ,b 的m a t l a b 程序如下： function nihex=[1.9,2.0,2.1,2.5,2.7,2.7,3.5,3.5,4.0,4.0,4.5,4.6,5.0,5.2,6.0,6.3,6.5,7.1,8.0,8.0,8.9,9.0,9.5,10.0];y=[1.4,1.3,1.8,2.5,2.8,2.5,3.0,2.7,4.0,3.5,4.2,3.5,5.5,5.0,5.5,6.4,6.0,5.3,6.5,7.0,8.5,8.0,8.1,8.1]; n=length(x); X=0; XY=0; XP=0; Y=0;for i=1:nX=X+x(i); Y=Y+y(i);XP=XP+x(i).*x(i); XY=XY+x(i)*y(i); endA=[24,X;X,XP]; B=[Y,XY]'; Gauss1(A,B)高斯消去法function x=Gauss1(A,b) n=length(b);for k=1:n-1if A(k,k)==0warning('系数矩阵奇异！'); return; endfor i=k+1:nL(i,k)=A(i,k)/A(k,k);A(i,k:n)=A(i,k:n)-L(i,k)*A(k,k:n); b(i)=b(i)-L(i,k)*b(k); end endif A(n,n)==0warning('系数矩阵奇异！'); return; endfor k=n:-1:1 if k==nx(n)=b(k)/A(n,n); elsex(k)=(b(k)-sum(A(k,k+1:n).*x(k+1:n)))/A(k,k); end end病态方程组的特例： 例：⎩⎨⎧=+=+0001.20001.12b a b a ⎩⎨⎧=+=+20001.12b a b a 比较两方程的准确解，可以发现它们的准确解差别很大下面估计A ∆和b ∆很小时解的误差 x x x ~-=∆x~和x 分别满足方程组b Ax b b x x A A =∆-=∆-∆-,))((两式相减，可得Ax b x A A ∆-∆=∆∆-)(当 A ∆很小时A A ∆-1很小，)(1A A E A A A ∆-=∆-- 可逆，)()(1Ax b A A x ∆-∆∆-=∆-)()(111Ax b A A A I x ∆-∆∆-=∆---)(1)()(11111x A b AA A Ax b A A A I ∆+∆∆-≤∆+∆∆-≤-----注意⎪⎪⎭⎫ ⎝⎛∆+∆∆-≤⎪⎪⎭⎫⎝⎛∆+∆∆-≤∆≤=∆≤∆------A A b b AA AA A A A AA x A b A AA A x xx A Ax b A A A A 11111111,,便知令⎪⎪⎭⎫⎝⎛∆+∆∆-≤∆==-A A b b AA k kx x x A A A cond k 1~,)(1误差估计式则得近似解此时表明，A A∆很小时，解的相对误差约为A ~和b~相对误差的k 倍；k 很大时，即使A ~和b ~的相对误差很小，解的相对误差也可能很大。

数值分析课程实验报告 实验名称 病态线性方程组的算法设计
班级
学号 姓名 序号
任课教师 评分 一、 实验目的
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)设计流程图：。

1 / 8数值分析实验六：解线性方程组的迭代法2016113 张威震1 病态线性方程组的求解1.1 问题描述理论的分析表明，求解病态的线性方程组是困难的。

实际情况是否如此，会出现怎样的现象呢？实验内容：考虑方程组Hx=b 的求解，其中系数矩阵H 为Hilbert 矩阵，,,1(),,,1,2,,1i j n n i j H h h i j n i j ⨯===+-这是一个著名的病态问题。

通过首先给定解（例如取为各个分量均为1）再计算出右端b 的办法给出确定的问题。

实验要求：（1）选择问题的维数为6，分别用Gauss 消去法、列主元Gauss 消去法、J 迭代法、GS 迭代法和SOR 迭代法求解方程组，其各自的结果如何？将计算结果与问题的解比较，结论如何？（2）逐步增大问题的维数（至少到100），仍然用上述的方法来解它们，计算的结果如何？计算的结果说明了什么？（3）讨论病态问题求解的算法1.2 算法设计首先编写各种求解方法的函数，Gauss 消去法和列主元高斯消去法使用实验5中编写的函数myGauss.m 即可，Jacobi 迭代法函数文件为myJacobi.m ，GS 迭代法函数文件为myGS.m ，SOR 方法的函数文件为mySOR.m 。

1.3 实验结果1.3.1 不同迭代法球求解方程组的结果比较选择H 为6*6方阵，方程组的精确解为x* = (1, 1, 1, 1, 1, 1)T ，然后用矩阵乘法计算得到b ，再使用Gauss 顺序消去法、Gauss 列主元消去法、Jacobi 迭代法、G-S 迭代法和SOR 方法分别计算得到数值解x1、x2、x3、x4，并计算出各数值解与精确解之间的无穷范数。

Matlab 脚本文件为Experiment6_1.m 。

迭代法的初始解x 0 = (0, 0, 0, 0, 0, 0)T ，收敛准则为||x(k+1)-x(k)||∞<eps=1e-6，SOR方法的松弛因子选择为w=1.3，计算结果如表1。