数值计算方法实验报告

重庆交通大学学生实验报告实验课程名称数值计算方法I开课实验室数学实验室学院理学院年级11专业班信息与计算科学学生姓名李伟凯学号631122020203开课时间2013至2014学年第1学期评分细则评分报告表述的清晰程度和完整性（20分）程序设计的正确性（40分）实验结果的分析（30分）实验方法的创新性（10分）总成绩教师签名邹昌文实验五解线性方程组的直接方法实验5.1（主元的选取与算法的稳定性）问题提出：Gauss消去法是我们在线性代数中已经熟悉的。但由于计算机的数值运算是在一个有限的浮点数集合上进行的，如何才能确保Gauss消去法作为数值算法的稳定性呢？Gauss消去法从理论算法到数值算法，其关键是主元的选择。主元的选择从数学理论上看起来平凡，它却是数值分析中十分典型的问题。实验内容：考虑线性方程组nnnRbRAbAx,,编制一个能自动选取主元，又能手动选取主元的求解线性方程组的Gauss消去过程。实验要求：（1）取矩阵1415157,6816816816bA，则方程有解Tx)1,,1,1(*。取n=10计算矩阵的条件数。让程序自动选取主元，结果如何？（2）现选择程序中手动选取主元的功能。每步消去过程总选取按模最小或按模尽可能小的元素作为主元，观察并记录计算结果。若每步消去过程总选取按模最大的元素作为主元，结果又如何？分析实验的结果。（3）取矩阵阶数n=20或者更大，重复上述实验过程，观察记录并分析不同的问题及消去过程中选择不同的主元时计算结果的差异，说明主元素的选取在消去过程中的作用。（4）选取其他你感兴趣的问题或者随机生成矩阵，计算其条件数。重复上述实验，观察记录并分析实验结果。实验5.2（线性代数方程组的性态与条件数的估计）问题提出：理论上，线性代数方程组bAx的摄动满足bbAAAAAcondxx11)(矩阵的条件数确实是对矩阵病态性的刻画，但在实际应用中直接计算它显然不现实，因为计算1A通常要比求解方程bAx还困难。实验内容：MATLAB中提供有函数“condest”可以用来估计矩阵的条件数，它给出的是按1-范数的条件数。首先构造非奇异矩阵A和右端，使得方程是可以精确求解的。再人为地引进系数矩阵和右端的摄动bA和，使得bA和充分小。实验要求：（1）假设方程Ax=b的解为x，求解方程bbxAAˆ)(，以1-范数，给出xxxxxˆ的计算结果。（2）选择一系列维数递增的矩阵（可以是随机生成的），比较函数“condest”所需机器时间的差别.考虑若干逆是已知的矩阵，借助函数“eig”很容易给出cond2(A)的数值。将它与函数“cond(A,2)”所得到的结果进行比较。（3）利用“condest”给出矩阵A条件数的估计，针对（1）中的结果给出xx的理论估计，并将它与（1）给出的计算结果进行比较，分析所得结果。注意，如果给出了cond(A)和A的估计，马上就可以给出1A的估计。（4）估计著名的Hilbert矩阵的条件数。njijihhHjinnji,,2,1,,11,)(,,思考题一：（Vadermonde矩阵）设niinniiniiniinnnnnnnxxxxbxxxxxxxxxxxxA002010022222121102001111，，其中，nkkxk,,1,0,1.01，（1）对n=2,5,8，计算A的条件数；随n增大，矩阵性态如何变化？（2）对n=5，解方程组Ax=b；设A的最后一个元素有扰动10-4，再求解Ax=b（3）计算（2）扰动相对误差与解的相对偏差，分析它们与条件数的关系。（4）你能由此解释为什么不用插值函数存在定理直接求插值函数而要用拉格朗日或牛顿插值法的原因吗？相关MATLAB函数提示：zeros(m,n)生成m行，n列的零矩阵ones(m,n)生成m行，n列的元素全为1的矩阵eye(n)生成n阶单位矩阵rand(m,n)生成m行,n列(0,1)上均匀分布的随机矩阵diag(x)返回由向量x的元素构成的对角矩阵tril(A)提取矩阵A的下三角部分生成下三角矩阵triu(A)提取矩阵A的上三角部分生成上三角矩阵rank(A)返回矩阵A的秩det(A)返回方阵A的行列式inv(A)返回可逆方阵A的逆矩阵[V，D]=eig(A)返回方阵A的特征值和特征向量norm(A,p)矩阵或向量的p范数cond(A,p)矩阵的条件数[L，U，P]=lu(A)选列主元LU分解R=chol(X)平方根分解Hi=hilb(n)生成n阶Hilbert矩阵实验程序：M文件程序为：functionx=gauss(n,r)n=input('请输入矩阵A的阶数:n=')A=diag(6*ones(1,n))+diag(ones(1,n-1),1)+diag(8*ones(1,n-1),-1)b=A*ones(n,1)p=input('条件数对应的范数是p-范数：p=')pp=cond(A,p)pause[m,n]=size(A);nb=n+1;Ab=[Ab]r=input('请输入是否为手动，手动输入1，自动输入0：r=')fori=1:n-1ifr==0[pivot,p]=max(abs(Ab(i:n,i)));ip=p+i-1;ifip~=iAb([iip],:)=Ab([ipi],:);disp(Ab);pauseendendifr==1i=iip=input('输入i列所选元素所处的行数：ip=');Ab([iip],:)=Ab([ipi],:);disp(Ab);pauseendpivot=Ab(i,i);fork=i+1:nAb(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-1:-1:1x(i)=(Ab(i,nb)-Ab(i,i+1:n)*x(i+1:n))/Ab(i,i);end（1）⑴取矩阵A的阶数:n=10，自动选取主元：formatlonggauss请输入矩阵A的阶数:n=10n=10条件数对应的范数是p-范数：p=1p=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=1r=1ans=1111111111②选取绝对值最小的元素为主元：gauss请输入矩阵A的阶数:n=10n=10条件数对应的范数是p-范数：p=2p=2pp=1.727556024913903e+003请输入是否为手动，手动输入1，自动输入0：r=1r=1ans=1.000000000000001.000000000000001.000000000000001.000000000000001.000000000000001.000000000000000.999999999999991.000000000000010.999999999999981.00000000000003（2）取矩阵A的阶数:n=10，手动选取主元：①选取绝对值最大的元素为主元：gauss请输入矩阵A的阶数:n=10n=10条件数对应的范数是p-范数：p=2p=2pp=1.727556024913903e+003请输入是否为手动，手动输入1，自动输入0：r=1r=1ans=1111111111②选取绝对值最小的元素为主元：gauss请输入矩阵A的阶数:n=10n=10条件数对应的范数是p-范数：p=2p=2pp=1.727556024913903e+003请输入是否为手动，手动输入1，自动输入0：r=1r=1ans=1.000000000000001.000000000000001.000000000000001.000000000000001.000000000000001.000000000000000.999999999999991.000000000000010.999999999999981.00000000000003（3）取矩阵A的阶数:n=20，手动选取主元：①选取绝对值最大的元素为主元：gauss请输入矩阵A的阶数:n=20条件数对应的范数是p-范数：p=1p=1pp=2.621437500000000e+006ans=11111111111111111111②选取绝对值最小的元素为主元：gauss请输入矩阵A的阶数:n=20.n=20条件数对应的范数是p-范数：p=2p=2pp=1.789670565881683e+006请输入是否为手动，手动输入1，自动输入0：r=1r=1ans=1.000000000000001.000000000000001.000000000000001.000000000000001.000000000000001.000000000000001.000000000000010.999999999999971.000000000000060.999999999999891.000000000000230.999999999999551.000000000000900.999999999998211.000000000003520.999999999993181.000000000012730.999999999978171.00000000002910（4）该题目的M程序如下所示functionx=gauss(n,r)n=input('请输入矩阵A的阶数:n=')A=hilb(n)b=A*ones(n,1)p=input('条件数对应的范数是p-范数：p=')pp=cond(A,p)pause[m,n]=size(A);nb=n+1;Ab=[Ab]r=input('请输入是否为手动，手动输入1，自动输入0：r=')fori=1:n-1ifr==0[pivot,p]=max(abs(Ab(i:n,i)));ip=p+i-1;ifip~=iAb([iip],:)=Ab([ipi],:);disp(Ab);pauseendendifr==1i=iip=input('输入i列所选元素所处的行数：ip=');Ab([iip],:)=Ab([ipi],:);disp(Ab);pauseendpivot=Ab(i,i);fork=i+1:nAb(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-1:-1:1x(i)=(Ab(i,nb)-Ab(i,i+1:n)*x(i+1:n))/Ab(i,i);end○1gauss请输入矩阵A的阶数：n=7n=7A=6100000861000008610000086100000861000008610000086b=7151515151514pp=3.174999999999999e+002Ab=61000007861000015086100015008610015000861015000086115000008614请输入是否为手动，手动输入1，自动输入0：r=1r=1条件数对应的范数是p-范数：p=1p=1i=1ans=0.999999999998691.000000000043370.999999999642991.0000000