数值分析-实验报告.docx

数值分析上机实验报告学院：专业：班级：学号：学生姓名：指导教师：实验五解线性方程组的直接方法实验5.1（主元的选取与算法的稳定性）问题提出：Gauss消去法是我们在线性代数中已经熟悉的。但由于计算机的数值运算是在一个有限的浮点数集合上进行的，如何才能确保Gauss消去法作为数值算法的稳定性呢？Gauss消去法从理论算法到数值算法，其关键是主元的选择。主元的选择从数学理论上看起来平凡，它却是数值分析中十分典型的问题。实验内容：考虑线性方程组Ax=b9AwRm编制一个能自动选取主元，又能手动选取主元的求解线性方程组的Gauss消去过程。实验要求：6 178 6 115(1)取矩阵4 = . .,b =:8 6 1158 614，那么方程有解N*= (1,1,1)7。取n=10计算矩阵的条件数。让程序自动选取主元，结果如何？(2)现选择程序中手动选取主元的功能。每步消去过程总选取按模最小或按模尽可能小的元素作为主元，观察并记录计算结果。假设每步消去过程总选取按模最大的元素作为主元，结果又如何？分析实验的结果。(3)取矩阵阶数n=20或者更大，重复上述实验过程，观察记录并分析不同的问题及消去过程中选择不同的主元时计算结果的差异，说明主元素的选取在消去过程中的作用。(4)选取其他你感兴趣的问题或者随机生成矩阵，计算其条件数。重复上述实验，观察记录并分析实验结果。实验过程：程序：建立M文件：functionx=gauss(n,r)n=inputC请输入矩阵A的阶数:n=)A=diag(6*ones(l,n)+diag(ones(1,n-1),1)+diag(8*ones(1,n-1),-1)b=A*ones(n,l)p=inp11C条件数对应的范数是P-范数：P=')pp=cond(A,p)pausem,n=size（八）;nb=n÷l;Ab=Abr=input(请输入是否为手动，手动输入1,自动输入0：r三,)fori=n-lifr=0pivot,p=max(abs(Ab(i:n,i);ip=p+i-l;ifip-=iAb(iip)=Ab(ipi,:);disp(Ab);pauseendendifr=li=iip=input(输入i列所选元素所处的行数：ip=，);Ab(iip)=Ab(ipi,:);disp(Ab);pauseendpivot=Ab(i,i);fork=i÷knAb(k,i:nb)=Ab(k,i:nb)-(Ab(k,i)/pivot)*Ab(i,i:nb);enddisp(Ab);pauseendx=zeros(n,1);x(n)=Ab(n,nb)/Ab(n,n);fori=n-l:-l:1x(i)=(Ab(i,nb)-Ab(i,i+lrn)*x(i+kn)Ab(i,i);end数值实验结果及分析：取矩阵A的阶数:n=10,自动选取主元：»formatlong»gauss请输入矩阵A的阶数:n=10n=10条件数对应的范数是P-范数：p=lP=1pp=2.557500000000000e+003请输入是否为手动，手动输入1,自动输入0：r=0r=0取矩阵A的阶数:n=10,手动选取主元：选取绝对值最大的元素为主元：»gauss请输入矩阵A的阶数:n=10n=10条件数对应的范数是P-范数：p=2P=2PP=1.727556024913903e+003请输入是否为手动，手动输入1,自动输入0：r=l=1ans=111111111选取绝对值最小的元素为主元：»gauss请输入矩阵A的阶数:n=10n=10条件数对应的范数是P-范数：p=2P=2pp=1.727556024913903e+003请输入是否为手动，手动输入1,自动输入0：r=l=1ans =1.000000000000001.000000000000000.999999999999991.000000000000001.000000000000001.000000000000011.000000000000001.000000000000000.999999999999981.00000000000003取矩阵A的阶数:n=20,手动选取主元: 选取绝对值最大的元素为主元:» gauss请输入矩阵A的阶数:20 条件数对应的范数是P-范数：p=lP= 1ans= 11111111Illl选取绝对值最小的元素为主元:» gauss请输入矩阵A的阶数:n=20.n =20条件数对应的范数是P-范数：p=2P= 2pp =1.789670565881683e+006请输入是否为手动，手动输入1,自动输入0： r=l= 1ans1.000000000000001.000000000000001.000000000000060.999999999999891.000000000000900.999999999993181.000000000000001.000000000000001.000000000000230.999999999998211.000000000012731.00000000000000 1.000000000000001.000000000000010.999999999999970.999999999999551.000000000003520.999999999978171.00000000002910将M文件中的第三行:A=diag(6*ones(1,n)+diag(ones(1,n-1),1)+diag(8*ones(1,n-1),-1)改为：A=hilb(n)>>gauss请输入矩阵A的阶数:n=7n=7条件数对应的范数是P-范数：p=lP=1请输入是否为手动，手动输入1,自动输入0：r=lr=1ans=1.000000000000510.999999999972511.000000000313540.999999998641331.000000002688050.999999997541811.00000000084337>>gauss请输入矩阵A的阶数:n=7n=7条件数对应的范数是P-范数：p=2P=2PP=4.753673569067072e+008请输入是否为手动，手动输入1,自动输入0：r=lr=1ans=0.999999999998691.000000000043370.999999999642991.000000001211430.999999998030381.000000001528250.99999999954491该问题在主元选取与算出结果有着很大的关系，取绝对值大的元素作为主元比取绝对值小的元素作为主元时产生的结果比拟准确，即选取绝对值小的主元时结果产生了较大的误差，条件数越大产生的误差就越大。讨论：在gauss消去法解线性方程组时，主元的选择与算法的稳定性有密切的联系，选取绝对值大的元素作为主元比绝对值小的元素作为主元时对结果产生的误差较小。条件数越大对用gauss消去法解线性方程组时，对结果产生的误差就越大。实验总结：对用gauss消去法解线性方程组时，主元的选取与算法的稳定性有密切的联系，选取适当的主元有利于得出稳定的算法，在算法的过程中，选取绝对值较大的主元比选取绝对值较小的主元更有利于算法的稳定，选取绝对值最大的元素作为主元时，得出的结果相对较准确较稳定。条件数越小，对用这种方法得出的结果更准确。在算除法的过程中要尽量防止使用较小的数做为除数，以免发生结果数量级加大，使大数吃掉小数，产生舍入误差。实验5.2(线性代数方程组的性态与条件数的估计)问题提出：理论上，线性代数方程组AX=人的摄动满足Cond（八）11A|AZ?|hW-1MlWWJ矩阵的条件数确实是对矩阵病态性的刻画，但在实际应用中直接计算它显然不现实，因为计算MIl通常要比求解方程Ax=A还困难。实验内容：MATLAB中提供有函数“condest”可以用来估计矩阵的条件数，它给出的是按I-范数的条件数。首先构造非奇异矩阵A和右端，使得方程是可以精确求解的。再人为地引进系数矩阵和右端的摄动A和助，使得|例网|明|充分小。实验要求：(1)假设方程Ax=b的解为X,求解方程(A+A)£="27,以I-范数，给出烙I=(U的计算结果。WWl(2)选择一系列维数递增的矩阵(可以是随机生成的)，比拟函数“CondeSt”所需机器时间的差异.考虑假设干逆是的矩阵，借助函数“eig”很容易给出Cond2（八）的数值。将它与函数Iond(A,2)”所得到的结果进行比拟。(3)利用“condes/给出矩阵A条件数的估计，针对(1)中的结果给出坤的理论估W计，并将它与(1)给出的计算结果进行比拟，分析所得结果。注意，如果给出了Cond（八）和网的估计，马上就可以给出卜的估计。(4)估计著名的HiIbert矩阵的条件数。H=(¼j)nxw,hij=-,Lj=1,2,2+«/一】实验过程：程序：n=input(,pleaseinputn:n-)a=fix(lOO*rand(n)+1x=ones(n,l)b=a*xdata=rand(n)*0.00001datb=rand(n,1)*0.00001A=a+dataB=b+datbxx=geshow(A,B)%输入矩阵的阶数%随机生成一个矩阵a%假设知道方程组的解全为1%用矩阵a和以知解得出矩阵b%随即生成扰动矩阵data%随即生成扰动矩阵datbx0=norm(xx-x,1)norm(x,l)%解扰动后的解%得出惇=%m的理论结果WM保存为：fanshu.mfunctionx=geshow(A,B)%用高斯消去法解方程组m,n=size（八）;nb=n+l;AB=AB;fori=l:n-lpivot=AB(i,i);fork=i+l:nAB(k,i:nb)=AB(k,i:nb)-(AB(k,i)/pivot)*AB(i,i：nb);endendx=zeros(n,l);x(n)=AB(n,nb)AB(n,n);fori=n-l:-l:lx(i)=(AB(i,nb)-AB(i,i+kn)*x(i+kn)AB(i,i);end保存为：geshow.mfunctioncond2（八）%自定义求二阶条件数B=A'*A;Vl,Dl=eig(B);V2,D2=eig(B(-l);cond2A=sqrt(max(max(D1)*sqrt(max(max(D2)end保存为：cond2.mformatlongforn=10:10:100n=nA=fix(100*randn(n);CondestA=Condest（八）cond2（八）condA2=cond(A,2)pauseend保存为：shiyan52.m%n为矩阵的阶%随机生成矩阵A%用condest求条件数%用自定义的求条件数%用cond求条件数%运行一次暂停%输入矩阵的阶数%随机生成一个矩阵a%假设知道方程组的解全为1%用矩阵a和以知解得出矩阵b%随即生成扰动矩阵data%随即生成扰动矩阵datb%利用第一小问的geshow.m求出解阵%得出焯=Rd的理论结果n=input(,pleaseinputn:n-)a=fix(lOO*rand(n)+1;x=ones(n,l);b=a*x;data=rand(n)*0.00001;datb=rand(n,1)*0.00001;A=a+data;B=b+datb;xx=geshow(A,B);xO=norm(xx-x,1)norm(x,1)xOO=cond（八）(l-norm(inv（八）)*norm(xx-x)*(norm(xx-x)(norm（八）)+norm(datb)norm(B)%得出datx=abs(x-x)保存为：Sy5_2.m?胃的估计值%求两者之间的误差formatlongforn=4:11n=nHi=hilb(n);cond1Hi=cond(Hi,l)cond2Hi=cond(Hi,2)condinfHi=cond(Hi,i11f)pauseend%n为矩阵的阶数%生成Hilbert矩阵%求Hilbert矩阵得三种条件数数值实验结果及分析:(1)»fanshupleaseinputn:n=66a=142516881989329385489260144088501316235219292324010100737241437227701x=111111b=251410221157218187data=1.0e-005*datb=1.0e-005*0.335632943352170.275099821466210.04452752039203A=1.0e+002*0.320000064986810.930000023756510.850000054592420.480000097047240.92000(X)35801710.60000002215793B=1.0e+002*xx=0.999998307797201.000000225695551.000000193415550.999999093880730.999999968940211.00000066032794x=胃号的计算结果为:(2)NCondestAcond2AcondA210e+00232.8905456307542132.89054563075420203.470959631940668e+002306.050503865112835e+002e+002e+002403.549487892582470e+00261.3753756968344861.3753756968336550e+002601.082004656409367e+0041.704830815154781e+0031.704830815108527e÷003703.234679145192132e+0033.878481155980936e+0023.878481155978439e+00280e+002902.063634143407935e+003e÷002e+0021001.536592818758897e+003e+002e+002(3)»sy5_2pleaseinputn:n=8n=8x=1.095033343195828e-006x00=1.705456352162135e-005datx=1.595953017842553e-005给出对悴=t三3的估计是：1.705456352162135e-005NIHI坤=1三1N的理论结果是：1.095033343195828e-006WHl结果相差：1.595953017842553e-005(4)nCondlHicond2HiCondinfHi42.837499999999738e+004e+0042.837499999999739e+00459.436559999999364e+0054.766072502414135e+0059.436559999999336e+00562.907027900294878e+007e+0072.907027900294064e÷0077e+0084.753673565864586e+008e+00883.387279082022742e+0101.525757545841988e+0103.387279081949470e+01091.099650993366047e÷012e+0111.099650991701052e÷012103.535372424347474e÷0131.602528637652488e+0133.535372455375642e÷013111.230369955362001e+0155.223946340715823e+0141.230369938308720e+015讨论：线性代数方程组的性态与条件数有着很重要的关系，既矩阵的条件数是刻画矩阵性质的一个重要的依据，条件数越大，矩阵“病态”性越严重，在解线性代数方程组的过程中较容易产生比拟大的误差，那么在实际问题的操作过程中，我们必须要减少对条件数来求解，把条件数较大的矩阵化成条件数较小的矩阵来进行求解。实验总结：在本次实验中，使我们知道了矩阵条件数对线性代数方程组求解的影响，条件数越大，对最后解的影响的越大，hilben矩阵是一个很“病态”的矩阵，他的条件数随着阶数的增加而增大，每增加一阶，条件数就增大一个数量级，在求解的过程中要尽量防止hilbert矩阵思考题一：(Vadermonde矩阵)设iX E 匕=()"I=OE(=o”z-5其中，xk=1+0.1c,A:=0,l,n,(1)对n=2,5,8,计算A的条件数；随n增大，矩阵性态如何变化？(2)对n=5,解方程组Ax=b;设A的最后一个元素有扰动10汽再求解A=b(3)计算(2)扰动相对误差与解的相对偏差，分析它们与条件数的关系。(4)你能由此解释为什么不用插值函数存在定理直接求插值函数而要用拉格朗日或牛顿插值法的原因吗？实验六解线性方程组的迭代法实验6.1(病态的线性方程组的求解)问题提出：理论的分析说明，求解病态的线性方程组是困难的。实际情况是否如此，会出现怎样的现象呢？实验内容：考虑方程组Hx=b的求解，其中系数矩阵H为HiIbert矩阵，H=(M)”xn，hij=.1,i,j=1,2,n这是一个著名的病态问题。通过首先给定解(例如取为各个分量均为1)再计算出右端b的方法给出确定的问题。实验要求：(1)选择问题的维数为6,分别用GaUSS消去法、J迭代法、GS迭代法和SOR迭代法求解方程组，其各自的结果如何？将计算结果与问题的解比拟，结论如何？12)逐步增大问题的维数，仍然用上述的方法来解它们，计算的结果如何？计算的结果说明了什么？(3)讨论病态问题求解的算法第问:选择问题的维数为6,分别用GaUSS消去法、J迭代法、GS迭代法和SoR迭代法求解方程组，其各自的结果如何？将计算结果与问题的解比拟，结论如何？源程序：Gauss消去法functionx=gauss(n)formatshortdispC请输入当前你要计算的Hilbert矩阵的阶数:)n=input(u);disp(构造出的HiIbert矩阵为:)A=hilb(n)d=diag(diag（八）),)*eye(n,n);b=A*ones(n,l);m,n=size（八）;nl=n+l;dispC该线性方程组的增广矩阵为1)Ab=Abfori=kn-l|zhuyuan,p=max(abs(Ab(i:n,i);ip=p+i-l;ifip=iAb(iip,:)=Ab(ipi,:);endzhuyuan=Ab(i,i);fork=i+l:nAb(k,i:n1)=Ab(k,i:n1)-(Ab(k,i)/ZhUyUan)*Ab(i,i:n1);enddisp(Ab);pauseendx=zeros(n,1);x(n)=Ab(n,n1)Ab(n,n);fori=n-l:-l:lx(i)=(Ab(i,n1)-Ab(i,i+kn)*x(i+kn)Ab(i,i);endJ迭代法formatshortdispC请输入当前你要计算的Hilbert矩阵的阶数:)n=input();disp(构造出的HiIbert矩阵为:)A=hilb(n)d=diag(diag（八）),)*eye(n,n);l=-(tril（八）-d);u=-(triu（八）-d);b=A*ones(n,l);x=0;0;0;0;0;0;fori=l:1000x1=inv(d)*(l+u)*x+inv(d)*b;x=x1;enddispC采用J迭代法结果为:，)xlGS迭代法fonnatshortdispC请输入当前你要计算的Hilbert矩阵的阶数:)n=input(");disp(构造出的HiIbert矩阵为:)A=hilb(n)d=diag(diag（八）)')*eye(n,n);l=-(tril（八）-d);u=-(triu（八）-d);b=A*ones(n,l);x=0;0;0;0;0;0;fori=kl000x1=inv(d-l)*u*x+inv(d-l)*b;x=x1;enddispC采用GS迭代法结果为:，)xlSOR迭代法formatshortdispC请输入你当前要计算的Hilbert矩阵的阶数:)n=input(");disp(构造出的Hilbert矩阵为:)A=hilb(n)d=diag(diag（八）)')*eye(n,n);l=-(tril（八）-d);u=-(triu（八）-d);b=A*ones(n,l);disp(请输入松弛因子Wr)w=input();x=0;0;0;0;0;0;fori=l:1000x1=inv(d-w*l)*(l-w)*d+w*u)*x+inv(d-w*l)*w*b;x=x1;enddispC采用SOR迭代法结果为:，)xl实验结果：Gauss消去法请输入当前你要计算的Hilbert矩阵的阶数：6构造出的Hilbert矩阵为：A=1.00000.50000.33330.25000.20000.16670.50000.33330.25000.20000.16670.14290.33330.25000.20000.16670.14290.12500.25000.20000.16670.14290.12500.11110.20000.16670.14290.12500.11110.10000.16670.14290.12500.11110.10000.0909该线性方程组的增广矩阵为：Ab=1.00000.50000.33330.25000.20000.16670.50000.33330.25000.20000.16670.14290.33330.25000.20000.16670.14290.12500.25000.20000.16670.14290.12500.11110.20000.16670.14290.12500.11110.10000.16670.14290.12500.11110.10000.09091.00000.50000.33330.25000.20000.166700.08330.08330.07500.06670.05952.45001.59291.21790.99560.84560.73652.45000.367900.08330.08890.08330.07620.06940.401200.07500.08330.08040.07500.06940.383100.06670.07620.07500.07110.06670.355600.05950.06940.06940.06670.06310.32821.00000.50000.33330.25000.20000.16672.450000.08330.08890.08330.07620.06940.401200-0.0056-0.0083-0.0095-0.0099-0.0333000.00330.00540.00640.00690.0221000.00510.00830.01020.01110.0347000.00600.00990.01220.01350.04161.00000.50000.33330.25000.20000.16672.450000.08330.08890.08330.07620.06940.4012000.00600.00990.01220.01350.0416000-0.0002-0.0004-0.0006-0.0013000-0.0001-0.0003-0.0004-0.00090000.00090.00190.00270.00551.00000.50000.33330.25000.20000.16672.450000.08330.08890.08330.07620.06940.4012000.00600.00990.01220.01350.04160000.00090.00190.00270.00550000-0.0000-0.0000-0.00010000-0.0000-0.0001-0.00011.00000.50000.33330.25000.20000.16672.450000.08330.08890.08330.07620.06940.4012000.00600.00990.01220.01350.04160000.00090.00190.00270.00550000-0.0000-0.0001-0.000100000-0.0000-0.0000ans=1.00001.00001.00001.00001.00001.0000J迭代法请输入当前你要计算的Hilbert矩阵的阶数:6构造出的Hilbert矩阵为:1.00000.50000.33330.25000.20000.16670.50000.33330.25000.20000.16670.14290.33330.25000.20000.16670.14290.12500.25000.20000.16670.14290.12500.11110.20000.16670.14290.12500.11110.10000.16670.14290.12500.11110.10000.0909采用J迭代法结果为：xl=NaNNaNNaNNaNNaNNaNGS迭代法dispC采用GS迭代法结果为:)xl请输入当前你要计算的Hilbert矩阵的阶数:6构造出的Hilbert矩阵为：A=1.00000.50000.33330.25000.20000.16670.50000.33330.25000.20000.16670.14290.33330.25000.20000.16670.14290.12500.25000.200.16670.14290.12500.11110.20000.16670.14290.12500.11110.10000.16670.14290.12500.11110.10000.0909采用GS迭代法结果为:xl=0.99931.01310.95371.03741.02960.9662SOR迭代法请输入你当前要计算的Hilbert矩阵的阶数:6构造出的HiIbert矩阵为：1.00000.50000.33330.25000.20000.16670.50000.33330.25000.20000.16670.14290.33330.25000.20000.16670.14290.12500.25000.20000.16670.14290.12500.11110.20000.16670.14290.12500.11110.10000.16670.14290.12500.11110.10000.0909请输入松弛因子w:1.5采用SOR迭代法结果为:xl=0.99851.02510.90761.10480.99180.9713结果分析整理以上采用四种方法求出的希尔伯特矩阵的解，如下表所示:GUESS迭代法J迭代法GS迭代法SOR迭代法1NAN0.99930.99851NAN1.01311.02511NAN0.95370.90761NAN1.03741.10481NAN1.02960.99181NAN0.96620.9713从上表可以看到，在希尔伯特矩阵为六阶的时候。采用GUESS迭代法可以得到精确解。在J迭代法时候迭代法发散，不能得出解，改良GUESS的GS迭代法得到了近似解，在采用松弛因子为1.5时，用SOR迭代法也得到了近似解。但是我们必不能就此得出GUESS迭代法最好等结论。为了得到更加合理的结果，接下来继续对不同阶矩阵进行讨论。第(2)问：逐步增大问题的维数，仍然用上述的方法来解它们，计算的结果如何？计算的结果说明了什么？源程序：GUESS迭代法依次改变高斯迭代法中矩阵的阶数，分别对N等于7到12运行程序，得到解。J迭代法formatshortforn=7:12A=hilb(n);d=diag(diag（八）),)*eye(n,n);l=-(tril（八）-d);u=-(triu（八）-d);b=A*ones(n,l);w=1.2;x=zeros(n,l);fori=l:l000x1=inv(d)*(l+u)*x+inv(d)*b;x=x1;endndispC采用J迭代法结果为:，)xlendGS迭代法formatshortforn=7:12A=hilb(n);d=diag(diag（八）),)*eye(n,n);l=-(tril（八）-d);u=-(triu（八）-d);b=A*ones(n,l);w=1.2;x=zeros(n,l);fori=1:1000x1=inv(d-l)*u*x+inv(d-l)*b;x=x1;endndispC采用GS迭代法结果为:)xlendSOR迭代法formatshortforn=7:12A=hilb(n);d=diag(diag（八）),)*eye(n,n);l=-(tril（八）-d);u=-(triu（八）-d);b=A*ones(n,l);w=1.2;x=zeros(n,l);fori=l:1000x1=inv(d-w*l)*(1-w)*d+w*u)*x+inv(d-w*l)*w*b;x=x1;endndispC采用SOR迭代法结果为xlend实验结果GUESS迭代法采用GUESS迭代法，分别计算阶数为7到12阶时，对应的解如下:7891011121.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00000.99991.00001.00001.0000