数值分析-实验报告.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.mformatlongfor
