数值实验 病态的线性方程组的求解
数值分析实验报告
实验五 解线性方程组的直接方法实验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矩阵)病态线性方程组的求解matlab程序
(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次作业(题目) 病态线性方程组的计算
第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法收敛快得多。
2.3方程组的病态问题
数值分析实验课第二单元第二部分数据拟合的数值方法实验报告问题重述:发现并观察数据拟合法正规方程组的病态问题.用拟合法求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)用 LU 分解和列主元高斯消去法求解上述方程组(③与⑥中的数据取小数 点后三位)输出 Ax=b 中矩阵 A 及向量 b,A=LU 分解的 L,U,detA 及解向量 x. (2)将方程组①中系数 3.01 改为 3.00,0.987 改为 0.990.用列主元高斯消 去法求解,输出列主元行交换次序、解向量 x 及 detA,并与(1)中结果比较. 将方程组②中的 2.099999 改为 2.1,5.900001 改为 5.9.用列主元高斯消去法求 解,输出解向量 x 及 detA,并与(1)中结果比较. 将③与⑥中的数据改为取小数点后四位用列主元高斯消去法求解, 输出解向 量 x 及 detA,并与(1)中结果比较. 判断④是否病态方程. 将⑤中的-4.273 改为-2.475 用列主元高斯消去法求解,输出解向量 x 及 detA,并与(1)中结果比较.
-2.000000 1.072000 5.643000 0 3.176000 1.801500 0 0 1.868071 det(A)为: -11.865990960000001 解向量 x 为: LU 分解解 x 为: x = -0.490396463271800 -0.051035181304353 0.367520253023993 (2) 列主元高斯消去法解 x 为: X = -0.490396463271872 -0.051035181304402 0.367520253024026 原方程是病态方程。
二、方法步骤
1.在 A 的非主对角线元素中,找出按模最大的元素 a pq ; 2.计算平面旋转矩阵 U pq ,其中的 sin( ) 及 cos( ) 由 cot 2
T AU pq , U U1U pq ( U1 的初始值取单位阵); 3.计算 A1 U pq (1) (1) ( 其中 aij 4. 如果 max aij 为 A1 的元素 ), 则停止计算 , 所求特征值为:
数值分析实验-病态线性方程组的算法设计
数值分析课程实验报告 实验名称 病态线性方程组的算法设计
班级
学号 姓名 序号
任课教师 评分 一、 实验目的
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)设计流程图:。
数值分析实验报告--实验6--解线性方程组的迭代法
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。
病态线性方程组的求解
《科学与工程计算》实验报告学号:姓名:一、实验内容:考虑方程组Hx=b 的求解,其中系数矩阵H 为Hilbert 矩阵,nj i j i h h H j i n n j i ,,2,1,,11,)(,, =-+==⨯这是一个著名的病态问题。
通过首先给定解(例如取为各个分量均为1)再计算出右端b 的办法给出确定的问题。
实验要求:(1)选择问题的维数为6,分别用Jacobi 迭代法、GS 迭代法和SOR 迭代法求解方程组,其各自的结果如何?将计算结果与问题的解比较,结论如何?(2)逐步增大问题的维数,仍然用上述的方法来解它们,计算的结果如何?计算的结果说明了什么? (3)讨论病态问题求解的算法。
二、程序设计的基本思想、原理和算法描述:1、 算法 Jacobi 迭代法若A 为稀疏矩阵,只需遍历非零元素GS 迭代法若A 为稀疏矩阵,只需遍历非零元素每步迭代计算量相当于一次矩阵与向量的乘法;不需要保留上一步的迭代解,与Jacobi 迭代法计算量一样。
SOR 迭代法(稠密矩阵)2、 函数组成double max(double array[100]) 求数组中的最大值函数3、 输入/输出设计对于方程组Hx=b 的求解,系数矩阵H 为Hilbert 矩阵,矩阵中的数由下列函数生成。
nj i j i h h H j i n n j i ,,2,1,,11,)(,, =-+==⨯X*取一个特解[1,1,1, (1)b 数组由矩阵的每行元素相加得来。
4、符号名说明double c[100] 用来存储第k+1次和第k 次迭代结果的差值的绝对值 double x[100] 第k+1次迭代结果,储存解数组 double x0[100] 初始向量double r 第k+1次和第k 次迭代结果的差值的绝对值的最大值 double sum 矩阵方程变换后右侧值的和 int k 迭代次数double a[100][100] 存储Hilbert 矩阵 double b[100] 存储b 向量三、源程序及注释:Jacobi 迭代法#include<iostream> #include<math.h> #include <iomanip> #include <stdio.h> using namespace std; int n;double max(double array[100])//求最大值函数 {double a=array[0]; int i;for(i=1; i<n; i++) {if(a<array[i]) a=array[i]; return a; } }double c[100]= {0.0};double x[100]= {0.0}; //第k+1次迭代结果,储存解数组 double x0[100]= {0.0}; //初始向量double r,sum=0; int main(){double s=0;int i,k,j;//k为迭代次数double a[100][100];double b[100]= {0.0};cout<<"请输入维数:"<<endl;cin>>n;cout<<"输出a数组:"<<endl;double max(double array[100]);for(i=0; i<n; i++){for(j=0; j<n; j++){a[i][j]=1.0/(i+j+1);printf("%10.6lf ",a[i][j]);//矩阵中的数精确到六位}cout<<endl;}cout<<"输出b数组:"<<endl;for(i=0; i<n; i++){for(j=0; j<n; j++){b[i]+=a[i][j];//矩阵每行的和}cout<<b[i]<<" ";}cout<<endl;cout<<"输入精度:"<<endl;cin>>s;for(k=1;; k++){for(i=0; i<n; i++){for(j=0; j<n; j++){sum=a[i][j]*x0[j]+sum;}x[i]=x0[i]+((b[i]-sum)/a[i][i]);//雅克比迭代方法计算方式c[i]=fabs(x[i]-x0[i]);//求差值的绝对值x0[i]=x[i];sum=0;r=max(c);if(r<s)//输出迭代次数{for(i=0; i<n; i++)cout<<"x"<<i<<" = "<<x[i]<<setprecision(10)<<endl;cout<<"迭代次数:"<<k<<endl;return 0;}}return 0;}GS迭代法#include<iostream>#include<math.h>#include <iomanip>#include <stdio.h>using namespace std;int n;double max(double array[100])//求最大值函数{double a=array[0];int i;for(i=1;i<n;i++){if(a<array[i])a=array[i];}return a;}main() {double s=0;double max(double array[100]);double c[100]={0.0};double x[100]={0.0};//第k+1次迭代结果,储存解数组double x0[100]={0.0};//初始向量int i,k,j;double r,sum=0;double a[100][100];double b[100]= {0.0};cout<<"请输入维数:"<<endl;cout<<"输出a数组:"<<endl;for(i=0; i<n; i++){for(j=0; j<n; j++){a[i][j]=1.0/(i+j+1);printf("%10.6lf ",a[i][j]);//矩阵中的数精确到六位}cout<<endl;}cout<<"输出b数组:"<<endl;for(i=0; i<n; i++){for(j=0; j<n; j++){b[i]+=a[i][j];//矩阵每行的和}cout<<b[i]<<" ";}cout<<endl;cout<<"输入精度:"<<endl;cin>>s;for(k=1;;k++){for(i=0;i<n;i++){for(j=0;j<n;j++){sum=a[i][j]*x0[j]+sum;}x[i]=x0[i]+((b[i]-sum)/a[i][i]);//gs迭代方法计算方式c[i]=fabs(x[i]-x0[i]);//求差值的绝对值x0[i]=x[i];sum=0;}r=max(c);if(r<s)//输出迭代次数{for(i=0;i<n;i++)cout<<"x"<<i<<" = "<<x[i]<<setprecision(10)<<endl;cout<<"迭代次数:"<<k<<endl;return 0;}}}SOR迭代法#include<iostream>#include<math.h>#include <iomanip>#include <stdio.h>using namespace std;int n;double max(double array[100]){double a=array[0];int i;for(i=1;i<n;i++){if(a<array[i])a=array[i];}return a;}int main() {double s=0,w=0;int i,k,j;//k为迭代次数double max(double array[100]);double c[100]= {0.0}; //double x[100]= {0.0}; //第k+1次迭代结果,储存解数组double x0[100]= {0.0}; //初始向量int i,k,j;double r,sum=0;double a[100][100];double b[100]= {0.0};cout<<"请输入维数:"<<endl;cin>>n;cout<<"输出a数组:"<<endl;for(i=0; i<n; i++){for(j=0; j<n; j++){a[i][j]=1.0/(i+j+1);printf("%10.6lf ",a[i][j]);//cout<<a[i][j]<<" ";}cout<<endl;}cout<<"输出b数组:"<<endl;for(i=0; i<n; i++){for(j=0; j<n; j++){b[i]+=a[i][j];}cout<<b[i]<<" ";}cout<<endl;cout<<"输入精度:"<<endl;cin>>s;cout<<"输入松弛因子:"<<endl;cin>>w;for(k=1;;k++){for(i=0;i<n;i++){for(j=0;j<n;j++){sum=a[i][j]*x0[j]+sum;}x[i]=x0[i]+(w*(b[i]-sum)/a[i][i]);//sor迭代方法的计算公式c[i]=fabs(x[i]-x0[i]);//求差值的绝对值x0[i]=x[i];sum=0;}r=max(c);if(r<s)//输出迭代次数{for(i=0;i<n;i++)cout<<"x"<<i<<" = "<<x[i]<<setprecision(10)<<endl; cout<<"迭代次数:"<<k<<endl;return 0;}}}四、运行输出结果:JacobiGSSOR五、结果比较分析:说明:Hx=b,H矩阵可以由函数hi,j=1/(i+j-1)直接由程序生成,为了设定参考解,我们先设x为分量全1的向量,求出b,然后将H和b作为已知量,求x,与设定的参考解对比。
病态线性方程组的求解
病态线性方程组的求解理论分析表明,数值求解病态线性方程组很困难。
考虑求解如下的线性方程组的求解Hx = b ,期中H 是Hilbert 矩阵,()ij n n H h ⨯=,11ij h i j =+-,i ,j = 1,2,…,n1. 估计矩阵的2条件数和阶数的关系2. 对不同的n ,取(1,1,,1)n x =∈ ,分别用Gauss 消去,Jacobi 迭代,Gauss-seidel 迭代,SOR 迭代和共轭梯度法求解,比较结果。
3. 结合计算结果,试讨论病态线性方程组的求解。
1)估计矩阵的2-条件数和阶数的关系矩阵的2-条件数定义为:1222()Cond A A A-=⨯,将Hilbert 矩阵带入有:1222()Cond H H H -=⨯。
使用MA TLAB 自带的cond2函数进行计算并画出log10(cond2)和阶数n 的关系曲线如下:可见当n 小于13的时候,条件数的对数与阶数有较好的线性关系,但是随着阶数的提高,条件数趋于“稳定”地振荡。
但是事实上,n较大时,H矩阵已经奇异,直接使用cond函数计算结果可能存在不准确性。
原因是对于条件数过大的矩阵使用inv函数求逆矩阵是不可靠的,应该使用invhilb函数进行求逆,并代入定义式中求解,生成的结果如下所示。
线性度较好,可知,Hilbert矩阵的2-条件数会随其阶数n的增加呈指数增大趋势,因此当n 较大时Hilbert矩阵将是严重病态的。
对不同的n,采用各种方法求解方程编写程序,选取n=2,3,5,10,15,20,迭代条件为迭代100000次或者是计算精度达到1e-6,若迭代次数少于设定最大值,表示相邻两次迭代达到精度要求,或者是计算结果超出范围。
X0取零向量,w取1.2,计算结果如下所示:由上可见,当n大于2时,Jacobi法已经不收敛,故其迭代次数已经没有意义。
当n=15时,Gauss消去法已经不收敛。
并且随着阶数的上升,gauss消去法的误差也随之上升。
数值求解Hilbert病态线性方程组
病态线性代数方程组的求解理论的分析表明,求解病态的线性代数方程组是困难的。
考虑方程组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 . 5 , 。 . 5 ) T , 于 是 = 可 I I x -  ̄ 1 = 。 5, 这表明解的相对误差是右端项相对误差的 1 。 0 0 倍.向
.
量 6 保 持 不 变 , 对 / 4 做 微 小 扰 动 ( _ 0 0 1 _ 0 0 . 0 . 0 0 1 1 ) ’ 计 算 得 _ 8 . 0 9 0 2 × 1 0 ~ , 则 扰 动 后 的 方 程 组
YONG L o ng -q ua n
( S c h o o l o f Ma t h e ma t i c s a n d C o m p u t e r S c i e n c e ,S h a a n x i U n i v e r s i t y o f T e c h n o l o g y ,H a n z h o n g 7 2 3 0 0 1 ,C h i n a)
9月
S e p . 2 0 1 5
文 章编 号 :1 0 0 7 — 9 8 3 1( 2 0 1 5)0 9 — 0 0 6 7 — 0 3
病态线性方程组算例分析
雍龙泉
( 陕西理工学院 数学与计算机科学学院 ,陕西 汉中 7 2 3 0 0 1 )
摘 要 :病 态线性 方程 组是 数值 分析 的 一个 重要研 究课 题 .给 出了 多个 病 态线性 方程组 的 算例 ,讨 论 了其病 态特性 ,介 绍 了求 解病 态线性 方程 组 的预 处理 方法 . 关键 词 :数值 分析 ;病 态线性 方程 组 ;条 件数 ;预 处理
Ab s t r a c t :S t u d y i n g i l l - c o n d i t i o n e d l i n e a r s y s t e m o f e q u a t i o n s i s a n i mp o r t a n t r e s e a r c h t o p i c i n n u me r i c a l
研究生数值分析9矩阵的条件数和病态线性方程组教学文案
Cond ( A)
即
Cond(A)A A1
通常使用的条件数有在行模意义下的条件数 与在谱模意义下的条件数,即
Cond(A)A1A
与
Cond(A)2A12
A 2
由线性代数知识,有
Cond(A)2
max(ATA) min(ATA)
定义 设A是非奇异矩阵,若 Co(nA)d 1
(1)仅b有小扰动δb
设方程组 AX=b+δb 的解为 X ~XX
即 A (XX)bb
②
②-①得 AXb
即 XA1b 于是有 X A1 b ③
另一方面,由①得
bA X 且 X 0
故
1A Xb
④
由③与④有
X A1
b
A⑤Leabharlann Xb表明解的相对误差不超过右端向量b 的相对误差的 A 1 A 倍。
(2)仅有小扰动δA(设 A+ δA 仍可逆)
因此称方程组
xx11
x2 2 1.0001x2
2
为病态方程组。
利用定义判断一个方程组是否病态,需要计算 矩阵的条件数,从而涉及计算逆矩阵,极不方便。
在实际计算中,常通过一些容易得到的信息来推断 方程组是否病态。 例如,当出现下列情况之一时,方程组很可能病态:
(1)用选主元消去法消元中出现小主元; (2)系数行列式的绝对值相对地很小; (3)系数矩阵元素间在数量级上相差很大且无 一定规律; (4)出现了相对地很大的解。
研究生数值分析9矩阵的条件数 和病态线性方程组
1 矩阵的条件数与线性方程组的性态 由于方程组AX=b系数矩阵A与右端向量b
的初始数据微小变化引起解的很大变化,这样 的方程组称为“病态”方程组。
病态矩阵线性方程组的数值解法
cond(A)z ̄II A Il 2 ll A Il 2=
1986。(01):1-10.
解 得 cond (A)z=8.0007 ×lo8 可证明矩阵 A病态程度非 常高 。
假记
=
:
(31)
. 32 ·
即 cond(A)),l是可判断矩 阵 A为病 态 可以是任一个秩 r的 m x n阶矩阵 ,且令 不 再举 例。
X-V,那么
参考文献
矩 阵 ;当 A的条件 数相对小 时 ,可判 断 A
【1】李庆扬 ,王能超 ,易大义。数 值分析【M1.4
是 良态矩 阵。 A的条件数越 大 ,其 方程 组
『I O 0]I
枷
1
...
^l l :
....
Seidel迭代
O
o
为
0
雾口 弓 )一赛 ))/乱 Ax ob,
显然 ,它 的精确解 为 【2 (1 ’ Fra bibliotek“》
一
1
程组解套 的影 响
b化 对方
即 , A(x+8 x)=b+8
(22) 因此,有 i ,
或
. ][薹]=[2 ·
版 。武 汉 :华 中科 技 大 学 出版 社 ,2006:
的病态程度越严重 。
VrATAV:A: O
0
185—212.
例 l’计算例 1中矩阵 A的谱条件数。
…
…
…
一
…
…
…
…
… · ^ … -
0
;0
[2】孔 祥元 ,魏 克让。线性 方程 组解的稳 定 性与病 态方程组的 最优 解法【J】.武测科技 ,
数学实验——线性代数方程组的数值解
实验5 线性代数方程组的数值解法分1 黄浩 2011011743一、实验目的1.学会用MATLAB软件数值求解线性代数方程组,对迭代法的收敛性和解的稳定性作初步分析;2.通过实例学习用线性代数方程组解决简化的实际问题。
二、实验内容1.《数学实验》第二版(问题1)问题叙述:通过求解线性方程组A1x=b1,A2x=b2,理解条件数的意义和方程组性态对解的影响,其中A1是n阶范德蒙矩阵,即A1=[1x0x02…x0n−11x1x12…x1n−1……………1x n−1x n−12 (x)n−1n−1], x k=1+0.1k , k=0,1,…,n−1A2是n阶希尔伯特矩阵,b1,b2分别是A1,A2的行和。
(1)编程构造A1(A2可直接用命令产生)和b1,b2;你能预先知道方程组A1x=b1和A2x=b2的解吗?令n=5,用左除命令求解(用预先知道的解可验证程序)。
(2)令n=5,7,9,…,计算A1和A2的条件数。
为观察他们是否病态,做以下试验:b1,b2不变,A1和A2的元素A1(n,n),A2(n,n)分别加扰动ε后求解;A1和A2不变,b1,b2的分量b1(n),b2(n)分别加扰动ε后求解。
分析A与b的微小扰动对解的影响。
ε取10^-10,10^-8,10^-6。
(3)经扰动得到的解记做x̃,计算误差‖x−x̃‖‖x‖,与用条件数估计的误差相比较。
模型转换及实验过程:(1)小题.由b1,b2为A1,A2的行和,可知方程组A1x=b1和A2x=b2的精确解均为n 行全1的列向量。
在n=5的情况下,用matlab编程(程序见四.1),构造A1,A2和b1,b2,使用高斯消去法得到的解x1,x2及其相对误差e1,e2(使用excel计算而得)为:由上表可见,当n=5时,所得的解都接近真值,误差在10^-12的量级左右。
(2)小题分别取n=5,7,9,11,13,15,计算A1和A2的条件数c1和c2,(程序见四.2),结果如下:由上表可见,二者的条件数都比较大,可能是病态的。
数值分析之病态线性方程组的解法
(
+
)
|| x || 1− || A || || A−1 || || A || || b || || A ||
|| A ||
数值分析
当方程组的系数矩阵A或右端项b受到 扰动ΔA, Δb时,引起的解的相对误差完全 由
来决定,它刻画了方程组的解对原始数据 的敏感程度。
数值分析
矩阵的条件数
定义:对非奇异矩阵A,称乘积||A|| ||A-1||为矩 阵A的条件数,记作
数值分析
三类扰动
线性方程组Ax=b的扰动可分为如下三类: 1. 系数矩阵A精确,常数项b有误差Δb; 2. b精确,A有误差ΔA; 3. A和b均有误差 ΔA和 Δb.
问题:求解Ax=b时,A和b的误差对x有何影响?
数值分析
扰动误差分析
1. 系数矩阵A精确,常数项b有误差Δb;
A(x + x) = b + b x = A−1b
=
2
1 3
1 n
1
n+1
1 n
1 n+1
1 2 n −1
cond (H2) = 27
cond (H3) 748
cond (H6) = 2.9 106
注:现在用Matlab数学软件可以很方便求矩阵的条件数!
数值分析
病态、良态线性方程组
定义:对线性方程组Ax=b,若cond(A)相对很大, 则称Ax=b是病态的线性方程组,也称A为病态矩阵; 若cond(A)相对很小,则称Ax=b是良态的线性方程 组,也称A为良态矩阵。
1.
;
2. A非奇异,k≠0,则cond(kA)=cond(A);
3. A非奇异对称矩阵,则cond(A)2=|λ1/λn|; 4. A是正交矩阵,则cond(A)2=1; 5. A可逆,R正交,则
数值分析希尔伯特病态线性方程组
病态线性代数方程组的求解理论的分析表明,求解病态的线性代数方程组是困难的。
考虑方程组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、下载文档前请自行甄别文档内容的完整性,平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益,请联系客服反馈,我们会尽快为您处理(人工客服工作时间:9:00-18:30)。
实验要求
1 )选择问题的维数为6,分别用Gauss消去法(即LU分解) J迭代法、GS迭代法和SOR迭代法求解方程组,其各自的 结果如何?将计算结果与问题的解比较,结论如何?
2)逐步增大问题的维数,仍然用上述的方法来解它们, 计算的结果如何?计算的结果说明了什么?
病态的线性方程组的求解
邹昌文Βιβλιοθήκη 问题的提出• 理论的分析表明,求解病态的线性方程 组是困难的,实际情况是否如此,会出 现怎样的现象呢?
实验内容
考虑方程组Hx b的求解,其中系数矩阵H为 Hilbert矩阵, 1 H (hi , j ) nn , hi , j , i, j 1,2,, n i j 1 通过首先给定解(例如取为各个分量均为 ) 1 再计算出右端b的办法给出确定的问题。