数值实验报告

1、重 庆 交 通 大 学学 生 实 验 报 告实验课程名称 数值计算方法II 开课实验室 数学实验室 学 院 理学院 年级 09 专业班 学 生 姓 名 学 号 开 课 时 间 2012 至 2013 学年第 2 学期评分细则评分报告表述的清晰程度和完整性（20分）程序设计的正确性（40分）实验结果的分析（30分）实验方法的创新性（10分）总成绩教师签名邹昌文实验五 解线性方程组的直接方法实验5.1 （主元的选取与算法的稳定性）问题提出：Gauss消去法是我们在线性代数中已经熟悉的。但由于计算机的数值运算是在一个有限的浮点数集合上进行的，如何才能确保Gauss消去法作为数值算法的稳定性呢？Gau

2、ss消去法从理论算法到数值算法，其关键是主元的选择。主元的选择从数学理论上看起来平凡，它却是数值分析中十分典型的问题。实验内容：考虑线性方程组 编制一个能自动选取主元，又能手动选取主元的求解线性方程组的Gauss消去过程。实验要求：（1）取矩阵，则方程有解。取n=10计算矩阵的条件数。让程序自动选取主元，结果如何？（2）现选择程序中手动选取主元的功能。每步消去过程总选取按模最小或按模尽可能小的元素作为主元，观察并记录计算结果。若每步消去过程总选取按模最大的元素作为主元，结果又如何？分析实验的结果。（3）取矩阵阶数n=20或者更大，重复上述实验过程，观察记录并分析不同的问题及消去过程中选择不同的

3、主元时计算结果的差异，说明主元素的选取在消去过程中的作用。（4）选取其他你感兴趣的问题或者随机生成矩阵，计算其条件数。重复上述实验，观察记录并分析实验结果。实验过程：实验代码：建立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，

4、自动输入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); pauseend

5、x=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实验结果及分析：取矩阵A的阶数:n=10，自动选取主元： format long gauss请输入矩阵A的阶数:n=10n = 10条件数对应的范数是p-范数：p=1p = 1pp = 2.557500000000000e+003请输入是否为手动，手动输入1，自动输入0：r=0r = 0取矩阵A的阶数:n=10，手动选取主元：选取绝对值最大的元素为主元： gauss请输入矩阵A的阶数:n=10n = 1

6、0条件数对应的范数是p-范数：p=2p = 2pp= 1.727556024913903e+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.727556024913903e+003请输入是否为手动，手动输入1，自动输入0：r=1r = 1ans = 1.00000000000000 1.00000000000000 1.00000000000000 1.00000000000000 1.0000

7、0000000000 1.00000000000000 0.99999999999999 1.00000000000001 0.99999999999998 1.00000000000003取矩阵A的阶数:n=20，手动选取主元： 选取绝对值最大的元素为主元： gauss请输入矩阵A的阶数:n=20条件数对应的范数是p-范数：p=1p = 1pp = 2.621437500000000e+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-范数：

8、p=2p = 2pp = 1.789670565881683e+006请输入是否为手动，手动输入1，自动输入0：r=1r = 1ans = 1.00000000000000 1.00000000000000 1.00000000000000 1.00000000000000 1.00000000000000 1.00000000000000 1.00000000000001 0.99999999999997 1.00000000000006 0.99999999999989 1.00000000000023 0.99999999999955 1.00000000000090 0.9999999

9、9999821 1.00000000000352 0.99999999999318 1.00000000001273 0.99999999997817 1.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=1p = 1pp = 9.851948872610030e+008请输入是否为手动，手动输入1，自动输入0：r=1r = 1ans = 1.00000000000051

10、0.99999999997251 1.00000000031354 0.99999999864133 1.00000000268805 0.99999999754181 1.00000000084337 gauss请输入矩阵A的阶数:n=7n = 7条件数对应的范数是p-范数：p=2p = 2pp = 4.753673569067072e+008请输入是否为手动，手动输入1，自动输入0：r=1r = 1ans = 0.99999999999869 1.00000000004337 0.99999999964299 1.00000000121143 0.99999999803038 1.0000

11、0000152825 0.99999999954491该问题在主元选取与算出结果有着很大的关系，取绝对值大的元素作为主元比取绝对值小的元素作为主元时产生的结果比较准确，即选取绝对值小的主元时结果产生了较大的误差，条件数越大产生的误差就越大。实验总结：通过本次实验,在gauss消去法解线性方程组时，主元的选择与算法的稳定性有密切的联系，选取绝对值大的元素作为主元比绝对值小的元素作为主元时对结果产生的误差较小。条件数越大对用gauss消去法解线性方程组时，对结果产生的误差就越大。主元的选取与算法的稳定性有密切的联系，选取适当的主元有利于得出稳定的算法，在算法的过程中，选取绝对值较大的主元比选取绝对

12、值较小的主元更有利于算法的稳定，选取绝对值最大的元素作为主元时，得出的结果相对较准确较稳定。条件数越小，对用这种方法得出的结果更准确。实验5.2（线性代数方程组的性态与条件数的估计）问题提出：理论上，线性代数方程组的摄动满足 矩阵的条件数确实是对矩阵病态性的刻画，但在实际应用中直接计算它显然不现实，因为计算通常要比求解方程还困难。实验内容：MATLAB中提供有函数“condest”可以用来估计矩阵的条件数，它给出的是按1-范数的条件数。首先构造非奇异矩阵A和右端，使得方程是可以精确求解的。再人为地引进系数矩阵和右端的摄动，使得充分小。实验要求：（1）假设方程Ax=b的解为x，求解方程，以1-范

13、数，给出的计算结果。（2）选择一系列维数递增的矩阵（可以是随机生成的），比较函数“condest”所需机器时间的差别.考虑若干逆是已知的矩阵，借助函数“eig”很容易给出cond2(A)的数值。将它与函数“cond(A,2)”所得到的结果进行比较。（3）利用“condest”给出矩阵A条件数的估计，针对（1）中的结果给出的理论估计，并将它与（1）给出的计算结果进行比较，分析所得结果。注意，如果给出了cond(A)和的估计，马上就可以给出的估计。（4）估计著名的Hilbert矩阵的条件数。实验过程：实验代码：n=input(please input n:n=) %输入矩阵的阶数a=fix(100

14、*rand(n)+1 %随机生成一个矩阵ax=ones(n,1) %假设知道方程组的解全为1b=a*x %用矩阵a和以知解得出矩阵bdata=rand(n)*0.00001 %随即生成扰动矩阵datadatb=rand(n,1)*0.00001 %随即生成扰动矩阵datbA=a+dataB=b+datbxx=geshow(A,B) %解扰动后的解x0=norm(xx-x,1)/norm(x,1) %得出的理论结果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

15、=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);endfunction cond2(A) %自定义求二阶条件数B=A*A;V1,D1=eig(B);V2,D2=eig(B(-1);cond2A=sqrt(max(max(D1)*sqrt(max(max(D2)endformat longfor n=10:10:100n=n %n为矩

16、阵的阶 A=fix(100*randn(n); %随机生成矩阵A condestA=condest(A) %用condest求条件数 cond2(A) %用自定义的求条件数 condA2=cond(A,2) %用cond求条件数 pause %运行一次暂停endn=input(please input n:n=) %输入矩阵的阶数a=fix(100*rand(n)+1; %随机生成一个矩阵ax=ones(n,1); %假设知道方程组的解全为1b=a*x; %用矩阵a和以知解得出矩阵bdata=rand(n)*0.00001; %随即生成扰动矩阵datadatb=rand(n,1)*0.0000

17、1; %随即生成扰动矩阵datbA=a+data;B=b+datb;xx=geshow(A,B); %利用第一小问的geshow.m求出解阵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) %求两者之间的误差format longfor n=4:11 n=n %n为矩阵的阶数Hi=hilb(n); %生成Hilbert矩阵cond1Hi=cond(Hi,1) %求Hilber

18、t矩阵得三种条件数 cond2Hi=cond(Hi,2) condinfHi=cond(Hi,inf) pauseend实验结果及分析： fanshuplease input n:n=6n = 6a = 14 25 16 88 19 89 32 93 85 48 92 60 14 40 88 50 13 16 23 52 19 29 2 32 40 10 100 7 37 24 14 3 72 27 70 1x = 1 1 1 1 1 1b = 251 410 221 157 218 187data = 1.0e-005 * 0.39690379186910 0.78196184196050

19、0.63712194084590 0.82064368228574 0.66093213223947 0.51488031898783 0.64986813059250 0.23756508204022 0.54592415509902 0.97047237460911 0.35801711338781 0.22157934638561 0.08500060621463 0.19573076378328 0.84805722441693 0.48692499554190 0.93819943010121 0.72500937095222 0.76880950325876 0.263213915

20、17561 0.80209765848011 0.81746853554695 0.48766697476487 0.06824661097009 0.96970170497170 0.71378506459614 0.66830641006672 0.64157116784600 0.09099035774397 0.96412426837254 0.71479723187621 0.97759973943565 0.67098263396985 0.30634935951390 0.67383411686207 0.20765658836866datb = 1.0e-005 * 0.161

21、11822555138 0.63822138259275 0.00022817289162 0.33563294335217 0.27509982146621 0.04452752039203A = 1.0e+002 * 0.14000003969038 0.25000007819618 0.16000006371219 0.88000008206437 0.19000006609321 0.89000005148803 0.32000006498681 0.93000002375651 0.85000005459242 0.48000009704724 0.92000003580171 0.

22、60000002215793 0.14000000850006 0.40000001957308 0.88000008480572 0.50000004869250 0.13000009381994 0.16000007250094 0.23000007688095 0.52000002632139 0.19000008020977 0.29000008174685 0.02000004876670 0.32000000682466 0.40000009697017 0.10000007137851 1.00000006683064 0.07000006415712 0.37000000909

23、904 0.24000009641243 0.14000007147972 0.03000009775997 0.72000006709826 0.27000003063494 0.70000006738341 0.01000002076566B = 1.0e+002 * 2.51000001611182 4.10000006382214 2.21000000002282 1.57000003356329 2.18000002750998 1.87000000445275xx = 0.99999830779720 1.00000022569555 1.00000019341555 0.9999

24、9909388073 0.99999996894021 1.00000066032794x0 = 6.181368174725440e-007的计算结果为：6.181368174725440e-007（2）NcondestAcond2AcondA2101.152530883943102e+00232.8905456307542132.89054563075420203.470959631940668e+00265.5412238417896665.54122384178720306.050503865112835e+0021.126539755706398e+0021.126539755706

25、322e+002403.549487892582470e+00261.3753756968344861.37537569683365506.855018184779408e+00281.1213899375359481.12138993753482601.082004656409367e+0041.704830815154781e+0031.704830815108527e+003703.234679145192132e+0033.878481155980936e+0023.878481155978439e+002808.318226153918658e+00286.2381429985251

26、386.23814299853018902.063634143407935e+0032.120696380331705e+0022.120696380331079e+0021001.536592818758897e+0031.559132035738491e+0021.559132035738373e+002 sy5_2please input n:n=8n = 8x0 = 1.095033343195828e-006x00 = 1.705456352162135e-005datx = 1.595953017842553e-005给出对的估计是：1.705456352162135e-005的理

27、论结果是： 1.095033343195828e-006结果相差： 1.595953017842553e-005(4)ncond1Hicond2HicondinfHi42.837499999999738e+0041.551373873892786e+0042.837499999999739e+00459.436559999999364e+0054.766072502414135e+0059.436559999999336e+00562.907027900294878e+0071.495105864009243e+0072.907027900294064e+00779.8519488971947

28、00e+0084.753673565864586e+0089.851948897198483e+00883.387279082022742e+0101.525757545841988e+0103.387279081949470e+01091.099650993366047e+0124.931544439891016e+0111.099650991701052e+012103.535372424347474e+0131.602528637652488e+0133.535372455375642e+013111.230369955362001e+0155.223946340715823e+0141

29、.230369938308720e+015实验总结：在本次实验中，使我们知道线性代数方程组的性态与条件数有着很重要的关系，既矩阵的条件数是刻画矩阵性质的一个重要的依据，条件数越大，矩阵“病态”性越严重，在解线性代数方程组的过程中较容易产生比较大的误差，则在实际问题的操作过程中，我们必须要减少对条件数来求解，把条件数较大的矩阵化成条件数较小的矩阵来进行求解。实验六 解线性方程组的迭代法实验6.1（病态的线性方程组的求解）问题提出：理论的分析表明，求解病态的线性方程组是困难的。实际情况是否如此，会出现怎样的现象呢？实验内容：考虑方程组Hx=b的求解，其中系数矩阵H为Hilbert矩阵， 这是一个著

30、名的病态问题。通过首先给定解（例如取为各个分量均为1）再计算出右端b的办法给出确定的问题。实验要求：（1）选择问题的维数为6，分别用Gauss消去法、J迭代法、GS迭代法和SOR迭代法求解方程组，其各自的结果如何？将计算结果与问题的解比较，结论如何？（2）逐步增大问题的维数，仍然用上述的方法来解它们，计算的结果如何？计算的结果说明了什么？（3）讨论病态问题求解的算法Gauss消去法实验代码：n=input(系数矩阵的阶数:);A=hilb(n);%构造Hilbert矩阵if a=0; for i=1:(n-1); for j=(i+1):n; x=A(j,i)/A(i,i); for k=1:

31、(n+1); A(j,k)=A(j,k)-x*A(i,k); end; end;end; y(n)=A(n,n+1)/A(n,n);for i=2:n; y(n-i+1)=A(n-i+1,n+1); for j=1:(i-1); y(n-i+1)=y(n-i+1)-A(n-i+1,n-j+1)*y(n-j+1); end; y(n-i+1)=y(n-i+1)/A(n-i+1,n-i+1);end;yend;%手动控制消元次序if a=1; for i=1:(n-1); A %显示每步消元的结果 m=input(请选取作为主消元行的行号); for l=1:(n+1); c=A(i,l); A(

32、i,l)=A(m,l); A(m,l)=c; end; for j=(i+1):n; x=A(j,i)/A(i,i); for k=1:(n+1); A(j,k)=A(j,k)-x*A(i,k); end; end;end;y(n)=A(n,n+1)/A(n,n);for i=2:n; y(n-i+1)=A(n-i+1,n+1); for j=1:(i-1); y(n-i+1)=y(n-i+1)-A(n-i+1,n-j+1)*y(n-j+1); end; y(n-i+1)=y(n-i+1)/A(n-i+1,n-i+1);end;yend;J迭代法实验代码：n=input(系数矩阵的阶数:);A

33、=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

34、; x2(i)=0;end;b2=A2*a02;for i=1:100000; 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

35、=1:100000; 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迭代方法迭代出现发散的不稳定现象，无法求解；用GS迭代方法迭代不发散，能求

36、得解，但收敛非常缓慢，当迭代次数取得相当大（100000次）时解仍在精确解附近波动；用SOR迭代方法迭代不发散，能求得解，当松弛因子去1.25左右时收敛较GS迭代快一些，但仍非常缓慢。2 选择问题的维数为20时：用Gauss消去法求得的解与精确解相差很大，相差10的量级。取初始向量为0，用J迭代方法迭代发散，无法求解；取初始向量为0，用GS迭代方法迭代不发散，能求得解，但收敛非常缓慢，迭代100000次后，算得的值与精确值1相差0.001的量级。取初始向量为0，用SOR迭代方法迭代不发散，能求得解，但收敛非常缓慢。从上面的结果可以看出当病态问题的阶数升高时作为直接法的Gauss消去法又能求解变

37、成不能求解。而GS和SOR迭代法在阶数升高时仍能求解。但在阶数较低时直接法能求得精确解而迭代发却总存在一定的误差。可见直接法与迭代法各有各的优势与不足。3 关于病态问题的求解，找到可逆的对角阵和使方程组化为 理论上最好选择对角阵满足：。实验七 非线性方程求根实验7.1（迭代法、初始值与收敛性）实验目的：初步认识非线性问题的迭代法与线性问题迭代法的差别，探讨迭代法及初始值与迭代收敛性的关系。问题提出：迭代法是求解非线性方程的基本思想方法，与线性方程的情况一样，其构造方法可以有多种多样，但关键是怎样才能使迭代收敛且有较快的收敛速度。实验内容：考虑一个简单的代数方程 针对上述方程，可以构造多种迭代法

38、，如 在实轴上取初始值x0，请分别用迭代（7.1）-（7.3）作实验，记录各算法的迭代过程。实验要求：（1）取定某个初始值，分别计算（7.1）-（7.3）迭代结果，它们的收敛性如何？重复选取不同的初始值，反复实验。请自选设计一种比较形象的记录方式（如利用MATLAB的图形功能），分析三种迭代法的收敛性与初值选取的关系。（2）对三个迭代法中的某个，取不同的初始值进行迭代，结果如何？试分析迭代法对不同的初值是否有差异？（3）线性方程组迭代法的收敛性是不依赖初始值选取的。比较线性与非线性问题迭代的差异，有何结论和问题。相关MATLAB函数提示：x=fzero(fun,x0) 返回一元函数fun的一个

39、零点，其中fun为函数句柄，x0为标量时，返回在x0附近的零点；x0为向量a,b时，返回函数在a,b中的零点x,f,h=fsolve(fun,x0) 返回一元或多元函数x0附近fun的一个零点，其中fun为函数句柄，x0为迭代初值；f返回fun在x的函数值，应该接近0； h返回值如果大于0，说明计算结果可靠，否则不可靠实验过程：程序：clearclcs=input(请输入要运行的方程，运行第几个输入几s=);clfif s=1 %决定坐标轴的范围和初始值a=-1.5;b=2.5; y00=0; x00=input(请输入第一个函数的初值：x00=);elseif s=2a=0.1;b=6.5; y00=0; x00=input(请输入第二个函数的初值：x00=);elseif s=3a=0;b=2; y00=0; x00=input(请输入第三个函数的初值：x00=);endx=linspace(a,b,80);y0=x; %计算直线y=xy1=zxy7f(x,s); %计算迭代函数y=f(