数值计算方法实验报告(共35页).doc

精选优质文档-倾情为你奉上重 庆 交 通 大 学学 生 实 验 报 告实验课程名称 数值计算方法I 开课实验室 数学实验室 学 院 理学院 年级11专业班 信息与计算科学 学 生 姓 名 李伟凯 学 号 3 开 课 时 间 2013 至 2014 学年第 1 学期评分细则评分报告表述的清晰程度和完整性（20分）程序设计的正确性（40分）实验结果的分析（30分）实验方法的创新性（10分）总成绩教师签名邹昌文实验五 解线性方程组的直接方法实验5.1 （主元的选取与算法的稳定性）问题提出：Gauss消去法是我们在线性代数中已经熟悉的。但由于计算机的数值运算是在一个有限的浮点数集合上进行的，如何才能确保Gauss消去法作为数值算法的稳定性呢？Gauss消去法从理论算法到数值算法，其关键是主元的选择。主元的选择从数学理论上看起来平凡，它却是数值分析中十分典型的问题。实验内容：考虑线性方程组 编制一个能自动选取主元，又能手动选取主元的求解线性方程组的Gauss消去过程。实验要求：（1）取矩阵，则方程有解。取n=10计算矩阵的条件数。让程序自动选取主元，结果如何？（2）现选择程序中手动选取主元的功能。每步消去过程总选取按模最小或按模尽可能小的元素作为主元，观察并记录计算结果。若每步消去过程总选取按模最大的元素作为主元，结果又如何？分析实验的结果。（3）取矩阵阶数n=20或者更大，重复上述实验过程，观察记录并分析不同的问题及消去过程中选择不同的主元时计算结果的差异，说明主元素的选取在消去过程中的作用。（4）选取其他你感兴趣的问题或者随机生成矩阵，计算其条件数。重复上述实验，观察记录并分析实验结果。实验5.2（线性代数方程组的性态与条件数的估计）问题提出：理论上，线性代数方程组的摄动满足 矩阵的条件数确实是对矩阵病态性的刻画，但在实际应用中直接计算它显然不现实，因为计算通常要比求解方程还困难。实验内容：MATLAB中提供有函数“condest”可以用来估计矩阵的条件数，它给出的是按1-范数的条件数。首先构造非奇异矩阵A和右端，使得方程是可以精确求解的。再人为地引进系数矩阵和右端的摄动，使得充分小。实验要求：（1）假设方程Ax=b的解为x，求解方程，以1-范数，给出的计算结果。（2）选择一系列维数递增的矩阵（可以是随机生成的），比较函数“condest”所需机器时间的差别.考虑若干逆是已知的矩阵，借助函数“eig”很容易给出cond2(A)的数值。将它与函数“cond(A,2)”所得到的结果进行比较。（3）利用“condest”给出矩阵A条件数的估计，针对（1）中的结果给出的理论估计，并将它与（1）给出的计算结果进行比较，分析所得结果。注意，如果给出了cond(A)和的估计，马上就可以给出的估计。（4）估计著名的Hilbert矩阵的条件数。思考题一：（Vadermonde矩阵）设 ，其中，（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文件程序为：function x=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)pausem,n=size(A);nb=n+1;Ab=A br=input('请输入是否为手动，手动输入1，自动输入0：r=')for i=1:n-1 if r=0 pivot,p=max(abs(Ab(i:n,i); ip=p+i-1; if ip=i Ab(i ip,:)=Ab(ip i,:);disp(Ab); pause end end if r=1 i=i ip=input('输入i列所选元素所处的行数：ip='); Ab(i ip,:)=Ab(ip i,:);disp(Ab); pause end pivot=Ab(i,i); for k=i+1:n Ab(k,i:nb)=Ab(k,i:nb)-(Ab(k,i)/pivot)*Ab(i,i:nb); end disp(Ab); pauseendx=zeros(n,1);x(n)=Ab(n,nb)/Ab(n,n);for i=n-1:-1:1 x(i)=(Ab(i,nb)-Ab(i,i+1:n)*x(i+1:n)/Ab(i,i);end（1）取矩阵A的阶数:n=10，自动选取主元：>> format long>> gauss请输入矩阵A的阶数:n=10n = 10条件数对应的范数是p-范数：p=1p = 1pp = 2.0000e+003请输入是否为手动，手动输入1，自动输入0：r=0r = 0取矩阵A的阶数:n=10，手动选取主元：选取绝对值最大的元素为主元：>> gauss请输入矩阵A的阶数:n=10n = 10条件数对应的范数是p-范数：p=2p = 2pp= 1.3903e+003请输入是否为手动，手动输入1，自动输入0：r=1r = 1ans= 1 1 1 1 1 1 1 1 1 1选取绝对值最小的元素为主元：>> gauss请输入矩阵A的阶数:n=10n = 10条件数对应的范数是p-范数：p=2p = 2pp = 1.3903e+003请输入是否为手动，手动输入1，自动输入0：r=1r = 1ans = 1.000 1.000 1.000 1.000 1.000 1.000 0.999 1.001 0.998 1.003（2）取矩阵A的阶数:n=10，手动选取主元：选取绝对值最大的元素为主元：>> gauss请输入矩阵A的阶数:n=10n = 10条件数对应的范数是p-范数：p=2p = 2pp= 1.3903e+003请输入是否为手动，手动输入1，自动输入0：r=1r = 1ans= 1 1 1 1 1 1 1 1 1 1选取绝对值最小的元素为主元：>> gauss请输入矩阵A的阶数:n=10n = 10条件数对应的范数是p-范数：p=2p = 2pp = 1.3903e+003请输入是否为手动，手动输入1，自动输入0：r=1r = 1ans = 1.000 1.000 1.000 1.000 1.000 1.000 0.999 1.001 0.998 1.003（3）取矩阵A的阶数:n=20，手动选取主元： 选取绝对值最大的元素为主元：>> gauss请输入矩阵A的阶数:n=20条件数对应的范数是p-范数：p=1p = 1pp = 2.0000e+006ans = 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 选取绝对值最小的元素为主元：>> gauss请输入矩阵A的阶数:n=20.n = 20条件数对应的范数是p-范数：p=2p = 2pp = 1.1683e+006请输入是否为手动，手动输入1，自动输入0：r=1r = 1ans = 1.000 1.000 1.000 1.000 1.000 1.000 1.001 0.997 1.006 0.989 1.023 0.955 1.090 0.821 1.352 0.318 1.273 0.817 1.910（4）该题目的M程序如下所示function x=gauss(n,r)n=input('请输入矩阵A的阶数:n=')A= hilb(n)b=A*ones(n,1)p=input('条件数对应的范数是p-范数：p=')pp=cond(A,p)pausem,n=size(A);nb=n+1;Ab=A br=input('请输入是否为手动，手动输入1，自动输入0：r=')for i=1:n-1 if r=0 pivot,p=max(abs(Ab(i:n,i); ip=p+i-1; if ip=i Ab(i ip,:)=Ab(ip i,:);disp(Ab); pause end end if r=1 i=i ip=input('输入i列所选元素所处的行数：ip='); Ab(i ip,:)=Ab(ip i,:);disp(Ab); pause end pivot=Ab(i,i); for k=i+1:n Ab(k,i:nb)=Ab(k,i:nb)-(Ab(k,i)/pivot)*Ab(i,i:nb); end disp(Ab); pauseendx=zeros(n,1);x(n)=Ab(n,nb)/Ab(n,n);for i=n-1:-1:1 x(i)=(Ab(i,nb)-Ab(i,i+1:n)*x(i+1:n)/Ab(i,i);end >>gauss请输入矩阵A的阶数：n=7n = 7A = 6 1 0 0 0 0 0 8 6 1 0 0 0 0 0 8 6 1 0 0 0 0 0 8 6 1 0 0 0 0 0 8 6 1 0 0 0 0 0 8 6 1 0 0 0 0 0 8 6b = 7 15 15 15 15 15 14pp = 3.9999e+002Ab = 6 1 0 0 0 0 0 7 8 6 1 0 0 0 0 15 0 8 6 1 0 0 0 15 0 0 8 6 1 0 0 15 0 0 0 8 6 1 0 15 0 0 0 0 8 6 1 15 0 0 0 0 0 8 6 14请输入是否为手动，手动输入1，自动输入0：r=1r = 1条件数对应的范数是p-范数：p=1p = 1i = 1ans = 0.869 1.337 0.299 1.143 0.038 1.825 0.491显然的是，该问题在主元选取与算出结果有着很大的关系，取绝对值大的元素作为主元比取绝对值小的元素作为主元时产生的结果比较准确，即选取绝对值小的主元时结果产生了较大的误差，条件数越大产生的误差就越大实验体会：运用高斯消去法求解线性方程组问题的时候，主元的选取和相应的消去法的选取决定了该算法的稳定性，选取绝对值大的元素比选取绝对值比较小的元素作为主元时对结果产生的误差影响比较小。并且增加条件数反而对结果的误差产生更大的影响。并且在运算中要尽量避免出现运用小数作为除数，使数量级加大，令大数吃掉小数的情况发生。实验5.2（线性代数方程组的性态与条件数的估计）问题提出：理论上，线性代数方程组的摄动满足 矩阵的条件数确实是对矩阵病态性的刻画，但在实际应用中直接计算它显然不现实，因为计算通常要比求解方程还困难。实验内容：MATLAB中提供有函数“condest”可以用来估计矩阵的条件数，它给出的是按1-范数的条件数。首先构造非奇异矩阵A和右端，使得方程是可以精确求解的。再人为地引进系数矩阵和右端的摄动，使得充分小。实验要求：（1）假设方程Ax=b的解为x，求解方程，以1-范数，给出的计算结果。（2）选择一系列维数递增的矩阵（可以是随机生成的），比较函数“condest”所需机器时间的差别.考虑若干逆是已知的矩阵，借助函数“eig”很容易给出cond2(A)的数值。将它与函数“cond(A,2)”所得到的结果进行比较。（3）利用“condest”给出矩阵A条件数的估计，针对（1）中的结果给出的理论估计，并将它与（1）给出的计算结果进行比较，分析所得结果。注意，如果给出了cond(A)和的估计，马上就可以给出的估计。（4）估计著名的Hilbert矩阵的条件数。程序代码：1保存M文件名为：fanshu.mn=input('please input n:n=') a=fix(100*rand(n)+1 x=ones(n,1) b=a*x data=rand(n)*0.00001 datb=rand(n,1)*0.00001 A=a+dataB=b+datbxx=geshow(A,B) x0=norm(xx-x,1)/norm(x,1) 保存M文件名为：geshow.m function x=geshow(A,B) m,n=size(A);nb=n+1;AB=A B;for i=1:n-1 pivot=AB(i,i); for k=i+1:n AB(k,i:nb)=AB(k,i:nb)-(AB(k,i)/pivot)*AB(i,i:nb); endendx=zeros(n,1);x(n)=AB(n,nb)/AB(n,n);for i=n-1:-1:1 x(i)=(AB(i,nb)-AB(i,i+1:n)*x(i+1:n)/AB(i,i);end2保存M文件名为：cond2.mfunction cond2(A) B=A'*A;V1,D1=eig(B);V2,D2=eig(B(-1);cond2A=sqrt(max(max(D1)*sqrt(max(max(D2)end保存M文件为：shiyan52.mformat longfor n=10:10:100n=n A=fix(100*randn(n); condestA=condest(A) cond2(A) condA2=cond(A,2) pause end3保存M文件为：sy5_2.mn=input('please input n:n=') %输入矩阵的阶数a=fix(100*rand(n)+1; %随机生成一个矩阵ax=ones(n,1); %假设知道方程组的解全为1b=a*x; data=rand(n)*0.00001; datb=rand(n,1)*0.00001; A=a+data;B=b+datb;xx=geshow(A,B); x0=norm(xx-x,1)/norm(x,1) x00=cond(A)/(1-norm(inv(A)*norm(xx-x)*(norm(xx-x)/(norm(A)+norm(datb)/norm(B) datx=abs(x0-x00) （4）format longfor n=4:11 n=n Hi=hilb(n); cond1Hi=cond(Hi,1) cond2Hi=cond(Hi,2) condinfHi=cond(Hi,inf) pauseend实验结果及其分析：(1)>>fanshuplease input n:n=6n=6 a = 82 28 96 80 68 71 91 55 49 96 76 4 13 96 81 66 75 28 92 97 15 4 40 5 64 16 43 85 66 10 10 98 92 94 18 83x = 1 1 1 1 1 1b = 425 371 359 253 284 395data = 1.0e-005 * Columns 1 through 4 0.5817 0.9002 0.8073 0.8377 0.0861 0.7063 0.2361 0.2143 0.8355 0.4379 0.8578 0.6081 0.2909 0.8231 0.3675 0.6133 0.6398 0.0900 0.3841 0.9777 0.3009 0.1265 0.4631 0.1137 Columns 5 through 6 0.5653 0.3803 0.9269 0.8679 0.5142 0.9058 0.6686 0.3737 0.5799 0.3663 0.5445 0.1531datb = 1.0e-005 * 0.8817 0.4989 0.7228 0.4809 0.1208 0.6031A = Columns 1 through 4 82.6227 28.7881 96.8307 80.6814 91.4802 55.9009 49.6819 96.0522 13.0488 96.6045 81.0766 66.9586 92.0801 97.3961 15.6769 4.7266 64.3596 16.2008 43.0039 85.7504 10.4570 98.0106 92.7353 94.9393 Columns 5 through 6 68.0592 71.5294 76.1152 4.4429 75.0512 28.0054 40.7226 5.2541 66.2530 10.2560 18.4253 83.1790B = 1.0e+002 * 4.8483 3.2497 3.6362 2.8376 2.9525 3.8386xx = 1.1561 0.4732 1.3244 0.6705 0.9937 0.3164x0 =1.4392e-006的计算结果为：1.4392e-006(2)NcondestAcond2AcondA2101.3102e+00232.42132.420203.0668e+00265.96665.720306.2835e+0021.6398e+0021.6322e+002403.2470e+00261.44861.365506.9408e+00281.59481.482601.9367e+0041.4781e+0031.8527e+003703.2132e+0033.0936e+0023.8439e+002808.8658e+00286.51386.018902.7935e+0032.1705e+0022.1079e+0021001.8897e+0031.8491e+0021.8373e+002(3)sy5_2please input n:n=8n = 8x0 = 3.9900e-007x00 = 5.3296e-007datx = 1.3396e-007给出对的估计是5.3296e-007的理论结果是：3.9900e-007结果相差：1.3396e-007(4)ncond1Hicond2HicondinfHi42.9738e+0041.2786e+0042.9739e+00459.9364e+0054.4135e+0059.9336e+00562.4878e+0071.9243e+0072.4064e+00779.4700e+0084.4586e+0089.8483e+00883.2742e+0101.1988e+0103.9470e+01091.6047e+0124.1016e+0111.1052e+012103.7474e+0131.2488e+0133.5642e+013111.2001e+0155.5823e+0141.8720e+015讨论：线性代数方程组的性态与条件数有着很重要的关系，既矩阵的条件数是刻画矩阵性质的一个重要的依据，条件数越大，矩阵“病态”性越严重，在解线性代数方程组的过程中较容易产生比较大的误差，则在实际问题的操作过程中，我们必须要减少对条件数来求解，把条件数较大的矩阵化成条件数较小的矩阵来进行求解。实验体会：在本次实验中，使我们知道了矩阵条件数对线性代数方程组求解的影响，条件数越大，对最后解的影响的越大，hilbert矩阵是一个很”病态”的矩阵，他的条件数随着阶数的增加而增大，每增加一阶，条件数就增大一个数量级，在求解的过程中要尽量避免hilbert矩阵实验六 解线性方程组的迭代法实验6.1（病态的线性方程组的求解）问题提出：理论的分析表明，求解病态的线性方程组是困难的。实际情况是否如此，会出现怎样的现象呢？实验内容：考虑方程组Hx=b的求解，其中系数矩阵H为Hilbert矩阵， 这是一个著名的病态问题。通过首先给定解（例如取为各个分量均为1）再计算出右端b的办法给出确定的问题。实验要求：（1）选择问题的维数为6，分别用Gauss消去法、J迭代法、GS迭代法和SOR迭代法求解方程组，其各自的结果如何？将计算结果与问题的解比较，结论如何？（2）逐步增大问题的维数，仍然用上述的方法来解它们，计算的结果如何？计算的结果说明了什么？（3）讨论病态问题求解的算法程序代码：Gauss消去法程序：function x=gauss(n,r)n=input('请输入矩阵的阶数:n=')A=hilb(n);%构造Hilbert矩阵p=input('条件数对应的范数是p-范数：p=')pp=cond(A,p)pausem,n=size(A);nb=n+1;Ab=A br=input('请输入是否为手动，手动输入1，自动输入0：r=')for i=1:n-1 if r=0 pivot,p=max(abs(Ab(i:n,i); ip=p+i-1; if ip=i Ab(i ip,:)=Ab(ip i,:);disp(Ab); pause end end if r=1 i=i ip=input('输入i列所选元素所处的行数：ip='); Ab(i ip,:)=Ab(ip i,:);disp(Ab); pause end pivot=Ab(i,i); for k=i+1:n Ab(k,i:nb)=Ab(k,i:nb)-(Ab(k,i)/pivot)*Ab(i,i:nb); end disp(Ab); pauseendx=zeros(n,1);x(n)=Ab(n,nb)/Ab(n,n);for i=n-1:-1:1 x(i)=(Ab(i,nb)-Ab(i,i+1:n)*x(i+1:n)/Ab(i,i);endJ迭代法程序：n=input('系数矩阵的阶数:');A=hilb(n);%构造Hilbert矩阵for i=1:n; a0(i)=1; %给定解 x(i)=0;end;b=A*a0' %由给定的解算出相应的b%进行迭代for i=1:100; y=x; for j=1:n; x(j)=b(j)/A(j,j); for k=1:j-1; x(j)=x(j)-A(j,k)*y(k)/A(j,j);end;for k=j+1:n; x(j)=x(j)-A(j,k)*y(k)/A(j,j);end;end;end;xGS迭代程序：n=input('系数矩阵的阶数:');%对题中给定的矩阵进行消元A2=hilb(n);for i=1:n; a02(i)=1; x2(i)=0;end;b2=A2*a02'for i=1:; for j=1:n; x2(j)=b2(j)/A2(j,j); for k=1:j-1; x2(j)=x2(j)-A2(j,k)*x2(k)/A2(j,j);end;for k=j+1:n; x2(j)=x2(j)-A2(j,k)*x2(k)/A2(j,j);end;end;end;x2SOR迭代程序：n=input('系数矩阵的阶数:');ss=input('松弛因子：');%对题中给定的矩阵进行消元A3=hilb(n);for i=1:n; a03(i)=1; x3(i)=0;end;b3=A3*a03'for i=1:; for j=1:n; rc=x3(j); x3(j)=b3(j)/A3(j,j); for k=1:j-1; x3(j)=x3(j)-A3(j,k)*x3(k)/A3(j,j);end;for k=j+1:n; x3(j)=x3(j)-A3(j,k)*x3(k)/A3(j,j);end;x3(j)=(1-ss)*rc+ss*x3(j);end;end;x3实验结果及其分析：给定各分量为1的解，计算出右端作为问题。1 选择问题的维数为6时：用Gauss消去法求得的解与精确解一致；取初始向量为0，用J迭代方法迭代出现发散的不稳定现象，无法求解；取初始向量为0，用GS迭代方法迭代不发散，能求得解，但收敛非常缓慢，当迭代次数取得相当大（次）时解仍在精确解附近波动；取初始向量为0，用SOR迭代方法迭代不发散，能求得解，当松弛因子去1.25左右时收敛较GS迭代快一些，但仍非常缓慢。2 选择问题的维数为20时：用Gauss消去法求得的解与精确解相差很大，相差10的量级。取初始向量为0，用J迭代方法迭代发散，无法求解；取初始向量为0，用GS迭代方法迭代不发散，能求得解，但收敛非常缓慢，迭代次后，算得的值与精确值1相差0.001的量级。取初始向量为0，用SOR迭代方法迭代不发散，能求得解，但收敛非常缓慢。从上面的结果可以看出当病态问题的阶数升高时作为直接法的Gauss消去法又能求解变成不能求解。而GS和SOR迭代法在阶数升高时仍能求解。但在阶数较低时直接法能求得精确解而迭代发却总存在一定的误差。可见直接法与迭代法各有各的优势与不足。关于病态问题的求解，主要的方法是对原方程作某些预处理，以降低系数矩阵的条件数。可以采取将系数矩阵A的每一行本别乘上适当常数的方法。即找到可逆的对角阵和使方程组化为 理论上最好选择对角阵满足：实验七 非线性方程求根实验7.1（迭代法、初始值与收敛性）实验目的：初步认识非线性问题的迭代法与线性问题迭代法的差别，探讨迭代法及初始值与迭代收敛性的关系。问题提出：迭代法是求解非线性方程的基本思想方法，与线性方程的情况一样，其构造方法可以有多种多样，但关键是怎样才能使迭代收敛且有较快的收敛速度。实验内容：考虑一个简单的代数方程针对上述方程，可以构造多种迭代法，如 在实轴上取初始值x0，请分别用迭代（7.1）-（7.3）作实验，记录各算法的迭代过程。实验要求：（1）取定某个初始值，分别计算（7.1）-（7.3）迭代结果，它们的收敛性如何？重复选取不同的初始值，反复实验。请自选设计一种比较形象的记录方式（如利用Matlab的图形功能），分析三种迭代法的收敛性与初值选取的关系。（2）对三个迭代法中的某个，取不同的初始值进行迭代，结果如何？试分析迭代法对不同的初值是否有差异？（3）线性方程组迭代法的收敛性是不依赖初始值选取的。比较线性与非线性问题迭代的差异，有何结论和问题。实验过程：7.1程序：7.1.1 对于第一个迭代方程的程序：保存为：diedai71.hClccleara=-1.5;b=2.5; y00=0; x00=input('请输入第一个函数的初值：x00='