数值分析实验报告1

1、WORD格式数值分析实验报告实验一误差分析实验 1.1 病态问题实验目的：算法有“优与“劣之分，问题也有“好与“坏之别。对数值方法的研究而言， 所谓坏问题就是问题本身对扰动敏感者， 反之属于好问题。通过本实验可获得一个初步体会。数值分析的大局部研究课题中，如线性代数方程组、矩阵特征值问题、非线性方程及方程组等都存在病态的问题。 病态问题要通过研究和构造特殊的算法来解决，当然一般要付出一些代价 如耗用更多的机器时间、 占用更多的存储空间等。问题提出：考虑一个高次的代数多项式20p(x) ( x 1)( x 2) ( x 20)(x k )(1.1)k1显然该多项式的全部根为 1,2,20 共计

2、20 个，且每个根都是单重的。现考虑该多项式的一个扰动p( x)x190(1.2)其中是一个非常小的数。这相当于是对1.1中x19的系数作一个小的扰动。我们希望比较 1.1和 1.2根的差异，从而分析方程 1.1的解对扰动的敏感性。实验内容： 为了实现方便，我们先介绍两个Matlab 函数：“roots和“ poly。uroots(a)其中假设变量 a 存储 n+1 维的向量，那么该函数的输出 u 为一个 n 维的向量。设 a 的元素依次为 a1, a2 , , an 1，那么输出 u 的各分量是多项式方程a1 xna2 xn 1an xan 10的全部根；而函数bpoly(v)的输出 b 是

3、一个 n+1 维变量，它是以n 维变量 v 的各分量为根的多项式的系数。可见“ roots和“ poly是两个互逆的运算函数。ess0.000000001;vezeros(1,21);ve(2)ess;专业资料整理WORD格式0专业资料整理WORD格式数值分析实验报告roots ( poly (1: 20)ve)上述简单的 Matlab 程序便得到 1.2的全部根，程序中的“ ess即是 1.2中的 。实验要求：1选择充分小的 ess，反复进展上述实验，记录结果的变化并分析它们。如果扰动项的系数 很小，我们自然感觉 1.1和 1.2的解应当相差很小。计算中你有什么出乎意料的发现？说明有些解关于

4、如此的扰动敏感性如何？2将方程 1.2中的扰动项改成x18或其它形式，实验中又有怎样的现象出现？3选作局部请从理论上分析产生这一问题的根源。注意我们可以将方程 1.2写成展开的形式，p( x,)x20x190(1.3)同时将方程的解x 看成是系数的函数，考察方程的某个解关于的扰动是否敏感，与研究它关于 的导数的大小有何关系？为什么？你发现了什么现象，哪些根关于 的变化更敏感？思考题一：上述实验的改进在上述实验中我们会发现用roots 函数求解多项式方程的精度不高，为此你可以考虑用符号函数 solve 来提高解的准确度， 这需要用到将多项式转换为符号多项式的函数 poly2sym, 函数的具体使

5、用方法可参考 Matlab 的帮助。实验过程：程序：a=poly(1:20);rr=roots(a);for n=2:21nfor m=1:9ess=10(-6-m);ve=zeros(1,21);ve(n)=ess;r=roots(a+ve);-6-ms=max(abs(r-rr)endend利用符号函数：思考题一a=poly(1:20);专业资料整理WORD格式1专业资料整理WORD格式数值分析实验报告y=poly2sym(a);rr=solve(y)for n=2:21nfor m=1:8ess=10(-6-m);ve=zeros(1,21);ve(n)=ess;a=poly(1:20)

6、+ve;y=poly2sym(a);r=solve(y);-6-ms=max(abs(r-rr)endend数值实验结果及分析：format long-6-m-7-8-9-10n2 2.79722687478331 1.86753632021158 1.06052762380748 0.252731442190473 1.69376699767424 0.92310666706964 0.08471614569741 0.4080402640941140.854013934155360.199410220210610.03972935295834050.110311005388710.0429

7、65323628440060000700008000090000100000110000120000130000140000150000160000170000180000190000200000210000-6-m-11-12-13-14n专业资料整理WORD格式2专业资料整理WORD格式数值分析实验报告20.038776764393800.162565848682800.13322664013598030.021642583175460004000050000600007000080000900001000001100001200001300001400001500001600001700

8、00180000190000200000210000讨论：利用这种方法进展这类实验， 可以很准确的扰动敏感性的一般规律。 即当对扰动项的系数越来越小时， 对其多项式扰动的结果也就越来越小， 即扰动敏感性与扰动项的系数成正比， 扰动项的系数越大， 对其根的扰动敏感性就越明显， 当扰动的系数一定时， 扰动敏感性与扰动的项的幂数成正比， 扰动的项的幂数越高， 对其根的扰动敏感性就越明显。实验总结：利用 MATLAB 来进展病态问题的实验，虽然其得出的结果是有误差的，但是可以很容易的得出对一个屡次的代数多项式的其中某一项进展很小的扰动，对其多项式的根会有一定的扰动的，所以对于这类病态问题可以借助于 M

9、ATLAB 来进展问题的分析。学号： 06450210专业资料整理WORD格式3专业资料整理WORD格式数值分析实验报告*：万轩实验二插值法实验 2.1 多项式插值的振荡现象问题提出 ：考虑一个固定的区间上用插值逼近一个函数。显然拉格朗日插值中使用的节点越多，插值多项式的次数就越高。 我们自然关心插值多项式的次数增加时，L(x) 是否也更加靠近被逼近的函数。 龙格给出了一个极著名例子。 设区间-1 ， 1 上函数f(x)=1(1+25x2)实验内容： 考虑区间 -1 ，1 的一个等距划分，分点为：x(i)=-1+2i/n,i=0,1,2,n泽拉格朗日插值多项式为：L(x)=l(i)(x)/(1

10、+25x(j)2 ) i=0,1,n其中 l(i)(x), i=0,1,n,n 是 n 次拉格朗日插值基函数。实验要求：选择不断增大的分点数目 n=2,3 , ，画出 f(x) 及插值多项式函数 L(x) 在 -1 ，1 上的图象，比较分析实验结果。2选择其它的函数，例如定义在区间-5 ，5 上的函数h(x)=x/(1+x4) , g(x)=arctanx重复上述的实验看其结果如何。3区间 a,b 上切比雪夫点的定义为：xk=b+a/2+(b-a)/2)cos(2k-1)/(2(n+1),k=1,2,n+1以 x1,x2x(n+1) 为插值节点构造上述各函数的拉格朗日插值多项式，比较其结果。实

11、验过程：程序：多项式插值的震荡现象实验2.1for m=1:6subplot(2,3,m)%把窗口分割成 2*3 大小的窗口largrang(6*m)%对 largrang 函数进展运行if m=1title('longn=6')elseif m=2title('longn=12')elseif m=3title('longn=18')elseif m=4title('longn=24')elseif m=5title('longn=30')elseif m=6title('longn=36')专业

12、资料整理WORD格式4专业资料整理WORD格式数值分析实验报告end%对每个窗口分别写上标题为插值点的个数end保存为： chazhi.mfunction largrang(longn)mm=input('please input mm(运行第几个函数就输入mm 为几 ):mm=')if mm=1%d 表示定义域的边界值d=1;elseif mm=2|mm=3d=5;endx0=linspace(-d,d,longn); %x 的节点 if mm=1y0=1./(1.+25.*x0.2);elseif mm=2y0=x0./(1.+x0.4);elseif mm=3y0=ata

13、n(x0);endx=sym('x');n=length(x0); s=0.0;for k=1:np=1.0;for j=1:nif j=kp=p*(x-x0(j)/(x0(k)-x0(j);endends=p*y0(k)+s;endy=s;if mm=1ezplot('1/(1+25*x2)')elseif mm=2ezplot('x/(1+x4)')elseif mm=3ezplot('atan(x)')endhold onezplot(y,-d,d)hold off保存为： largrang.m数值实验结果及分析：对于第一个

14、函数f(x)=1/(1+25x2)专业资料整理WORD格式5专业资料整理WORD格式0.60.40.20-11.510.50-0.5-10.50-0.5-5420-2-4-5数 值分析实验报 告longn=6longn=12longn=18110.50.500-0.501-101-101xxxlongn=24longn=30longn=3621.51100.5-10-2-0.501-101-101xxx对于第二个函数h(x)=x/(1+x4)longn=6longn=12longn=181100-1-105-505-505xxxlongn=24longn=3020longn=36200-2-2

15、005-505-505xxx对于第三个函数 g(x)=arctan(x)专业资料整理WORD格式6专业资料整理WORD格式数 值分析实验报 告longn=6longn=12longn=181221000-1-1-2-2-505-505-505xxxlongn=24longn=30longn=36222000-2-2-2-505-505-505xxx讨论：通过对三个函数得出的 largrang 插值多项式并在数学软件中的运行， 得出函数图象，说明了对函数的支点不是越多越好，而是在函数的两端而言支点越多，而 largrang 插值多项式不是更加靠近被逼近的函数，反而更加远离函数，在函数两端的跳动性

16、更加明显， argrang 插值多项式对函数不收敛。实验总结：利用 MATLAB 来进展函数的 largrang 插值多项式问题的实验，虽然其得出的结果是有误差的，但是增加支点的个数进展屡次实验，可以找出函数的 largrang 插值多项式的一般规律，当支点增加时， largrang 插值多项式对函数两端不收敛，不是更加逼近，而是更加远离，跳动性更强。所以对于函数的 largrang 插值多项式问题可以借助于 MATLAB 来进展问题的分析，得到比较准确的实验结规律。学号： 06450210专业资料整理WORD格式7专业资料整理WORD格式数值分析实验报告*：万轩实验五解线性方程组的直接方法实

17、验 5.1 主元的选取与算法的稳定性问题提出： Gauss消去法是我们在线性代数中已经熟悉的。但由于计算机的数值运算是在一个有限的浮点数集合上进展的， 如何才能确保 Gauss消去法作为数值算法的稳定性呢？ Gauss消去法从理论算法到数值算法，其关键是主元的选择。主元的选择从数学理论上看起来平凡，它却是数值分析中十分典型的问题。实验内容：考虑线性方程组Ax b,ARn n , bRn编制一个能自动选取主元， 又能手动选取主元的求解线性方程组的Gauss消去过程。实验要求：61786115 1取矩阵A, b，那么方程有解 x*(1,1,1)T。取 n=10861158614计算矩阵的条件数。让

18、程序自动选取主元，结果如何？（ 2现选择程序中手动选取主元的功能。每步消去过程总选取按模最小或按模尽可能小的元素作为主元， 观察并记录计算结果。 假设每步消去过程总选取按模最大的元素作为主元，结果又如何？分析实验的结果。（ 3取矩阵阶数 n=20 或者更大，重复上述实验过程， 观察记录并分析不同的问题及消去过程中选择不同的主元时计算结果的差异， 说明主元素的选取在消去过程中的作用。（ 4选取其他你感兴趣的问题或者随机生成矩阵，计算其条件数。重复上述实验，观察记录并分析实验结果。实验过程：程序：建立 M 文件：function x=gauss(n,r)n=input('请输入矩阵 A 的

19、阶数 :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('条件数对应的X数是p-X数： p=')pp=cond(A,p)pausem,n=size(A);nb=n+1;Ab=A b专业资料整理WORD格式8专业资料整理WORD格式数值分析实验报告r=input(' 请输入是否为手动，手动输入1，自动输入 0： r=')for i=1:n-1if r=0pivot,p=max(abs(Ab(i:n,i);ip=p+i-1;if ip=

20、iAb(i ip,:)=Ab(ip i,:);disp(Ab); pauseendendif r=1i=iip=input(' 输入 i 列所选元素所处的行数：ip=');Ab(i ip,:)=Ab(ip i,:);disp(Ab); pauseendpivot=Ab(i,i);for k=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);for i=n-1:-1:1x(i)=(Ab(i,nb)-Ab(i,i+

21、1:n)*x(i+1:n)/Ab(i,i);end数值实验结果及分析：取矩阵 A 的阶数 :n=10，自动选取主元：>> format long>> gauss请输入矩阵 A 的阶数 :n=10n =10条件数对应的X数是p-X数： p=1p =1pp =2.557500000000000e+003请输入是否为手动，手动输入1，自动输入 0：r=0r =0取矩阵 A 的阶数 :n=10，手动选取主元：选取绝对值最大的元素为主元：>> gauss请输入矩阵 A 的阶数 :n=10n =10条件数对应的X数是p- X数： p=2p =2pp=1.72755602

22、4913903e+003专业资料整理WORD格式9专业资料整理WORD格式数值分析实验报告请输入是否为手动，手动输入 1，自动输入 0：r=1 r = 1ans=1111111111选取绝对值最小的元素为主元：>> gauss请输入矩阵 A 的阶数 :n=10n =10条件数对应的X数是p- X数： p=2p =2pp = 1.727556024913903e+003请输入是否为手动，手动输入 1，自动输入 0：r=1r =1ans =1.000000000000001.000000000000001.000000000000001.000000000000001.00000000

23、0000001.000000000000000.999999999999991.000000000000010.999999999999981.00000000000003取矩阵 A 的阶数 :n=20，手动选取主元： 选取绝对值最大的元素为主元：>> gauss请输入矩阵 A 的阶数 :n=20条件数对应的X数是p- X数： p=1p =1pp =2.621437500000000e+006ans =11111111111111111111 选取绝对值最小的元素为主元：>> gauss请输入矩阵 A 的阶数 :n=20.n =20条件数对应的X数是p- X数： p=2

24、p =2pp = 1.789670565881683e+006请输入是否为手动，手动输入 1，自动输入 0：r=1r =1ans =1.000000000000001.000000000000001.000000000000001.000000000000001.000000000000001.000000000000001.000000000000010.999999999999971.000000000000060.999999999999891.000000000000230.999999999999551.000000000000900.999999999998211.00000000

25、0003520.999999999993181.000000000012730.999999999978171.00000000002910将 M文件中的第三行：A=diag(6*ones(1,n)+diag(ones(1,n-1),1)+diag(8*ones(1,n-1),-1)专业资料整理WORD格式10专业资料整理WORD格式数值分析实验报告改为：A=hilb(n) >> gauss请输入矩阵 A 的阶数 :n=7n =7条件数对应的X数是p- X数： p=1p =1pp = 9.851948872610030e+008请输入是否为手动，手动输入 1，自动输入 0：r=1r

26、 =1ans =1.000000000000510.999999999972511.000000000313540.999999998641331.000000002688050.999999997541811.00000000084337 >> gauss请输入矩阵 A 的阶数 :n=7n =7条件数对应的X数是p- X数： p=2p =2pp = 4.753673569067072e+008请输入是否为手动，手动输入 1，自动输入 0：r=1r =1ans =0.999999999998691.000000000043370.999999999642991.0000000012

27、11430.999999998030381.000000001528250.99999999954491该问题在主元选取与算出结果有着很大的关系， 取绝对值大的元素作为主元比取绝对值小的元素作为主元时产生的结果比较准确， 即选取绝对值小的主元时结果产生了较大的误差，条件数越大产生的误差就越大。讨论：在 gauss消去法解线性方程组时， 主元的选择与算法的稳定性有密切的联系，选取绝对值大的元素作为主元比绝对值小的元素作为主元时对结果产生的误差较小。条件数越大对用 gauss消去法解线性方程组时， 对结果产生的误差就越大。实验总结：对用 gauss消去法解线性方程组时，主元的选取与算法的稳定性有密

28、切的联系，选取适当的主元有利于得出稳定的算法， 在算法的过程中， 选取绝对值较大的主元比选取绝对值较小的主元更有利于算法的稳定， 选取绝对值最大的元素作为主元时，得出的结果相对较准确较稳定。 条件数越小， 对用这种方法得出的结果更准确。在算除法的过程中要尽量防止使用较小的数做为除数， 以免发生结果数量级加大，使大数吃掉小数，产生舍入误差。学号： 06450210专业资料整理WORD格式11专业资料整理WORD格式数值分析实验报告*：万轩实验 5.2线性代数方程组的性态与条件数的估计问题提出：理论上，线性代数方程组Axb 的摄动满足xc o n(dA)Abx1 A 1AAb矩阵的条件数确实是对矩

29、阵病态性的刻画， 但在实际应用中直接计算它显然不现实，因为计算 A 1通常要比求解方程 Ax b 还困难。实验内容： Matlab 中提供有函数“ condest可以用来估计矩阵的条件数，它给出的是按 1-X数的条件数。首先构造非奇异矩阵 A 和右端，使得方程是可以精确求解的。再人为地引进系数矩阵和右端的摄动A和 b ，使得A 和 b 充分小。实验要求： 1假设方程 Ax=b 的解为 x，求解方程( A"A)x b b ，以1-X数，给出*"x的计算结果。x x（ 2选择一系列维数递增的矩阵可以是随机生成的 ，比较函数“ condest所需机器时间的差异 .考虑假设干逆是的

30、矩阵， 借助函数“eig很容易给出 cond2(A) 的数值。将它与函数“ cond(A,2)所得到的结果进展比较。 3利用“condest给出矩阵 A 条件数的估计， 针对 1中的结果给出x 的x理论估计，并将它与 1给出的计算结果进展比较，分析所得结果。注意，如果给出了 cond(A)和 A 的估计，马上就可以给出A1的估计。 4估计著名的 Hilbert 矩阵的条件数。H (hi , j ) n n , hi, j1, i , j 1,2, , nij 1实验过程：程序：n=input('please input n:n=')%输入矩阵的阶数a=fix(100*rand(

31、n)+1%随机生成一个矩阵ax=ones(n,1)%假设知道方程组的解全为1b=a*x%用矩阵 a 和以知解得出矩阵bdata=rand(n)*0.00001%随即生成扰动矩阵data专业资料整理WORD格式12专业资料整理WORD格式数值分析实验报告datb=rand(n,1)*0.00001%随即生成扰动矩阵 datbA=a+dataB=b+datb*=geshow(A,B)%解扰动后的解x0=norm(*-x,1)/norm(x,1)*"x%得出的理论结果xx保存为： fanshu.mfunction x=geshow(A,B)%用高斯消去法解方程组m,n=size(A);nb

32、=n+1;AB=A B;for i=1:n-1pivot=AB(i,i);for k=i+1:nAB(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:1x(i)=(AB(i,nb)-AB(i,i+1:n)*x(i+1:n)/AB(i,i);end保存为： geshow.mfunction cond2(A)%自定义求二阶条件数B=A'*A;V1,D1=eig(B);V2,D2=eig(B(-1);cond2A=sqrt(max(max(D

33、1)*sqrt(max(max(D2)end保存为： cond2.mformat longfor n=10:10:100n=n%n 为矩阵的阶A=fix(100*randn(n);%随机生成矩阵 AcondestA=condest(A)%用 condest求条件数cond2(A)%用自定义的求条件数condA2=cond(A,2)%用 cond 求条件数pause%运行一次暂停end保存为： shiyan52.m专业资料整理WORD格式13专业资料整理WORD格式数值分析实验报告n=input('please input n:n=')%输入矩阵的阶数a=fix(100*rand

34、(n)+1;%随机生成一个矩阵 ax=ones(n,1);%假设知道方程组的解全为 1b=a*x;%用矩阵 a 和以知解得出矩阵 bdata=rand(n)*0.00001;%随即生成扰动矩阵 datadatb=rand(n,1)*0.00001;%随即生成扰动矩阵 datbA=a+data;B=b+datb;*=geshow(A,B);%利用第一小问的 geshow.m 求出解阵x0=norm(*-x,1)/norm(x,1)*"x%得出的理论结果xxx00=cond(A)/(1-norm(inv(A)*norm(*-x)*(norm(*-x)/(norm(A)+norm(datb

35、)/norm(B)%得出x"x x 的估计值xxdatx=abs(x0-x00)%求两者之间的误差保存为： sy5_2.mformat longfor n=4:11n=n%n 为矩阵的阶数Hi=hilb(n);%生成 Hilbert 矩阵cond1Hi=cond(Hi,1)%求 Hilbert 矩阵得三种条件数cond2Hi=cond(Hi,2)condinfHi=cond(Hi,inf)pauseend数值实验结果及分析： >> fanshuplease input n:n=6n =6a =1425168819893293854892601440885013162352

36、19292324010100737241437227701x =专业资料整理WORD格式14专业资料整理WORD格式数值分析实验报告111111b =251410221157218187data =1.0e-005 *0.396903791869100.781961841960500.637121940845900.820643682285740.660932132239470.514880318987830.649868130592500.237565082040220.545924155099020.970472374609110.358017113387810.2215793463856

37、10.085000606214630.195730763783280.848057224416930.486924995541900.938199430101210.725009370952220.768809503258760.263213915175610.802097658480110.817468535546950.487666974764870.068246610970090.969701704971700.713785064596140.668306410066720.641571167846000.090990357743970.964124268372540.714797231

38、876210.977599739435650.670982633969850.30634935951390 0.67383411686207 0.20765658836866 datb =1.0e-005 *0.161118225551380.638221382592750.000228172891620.335632943352170.275099821466210.04452752039203A =1.0e+002 *0.140000039690380.250000078196180.160000063712190.880000082064370.190000066093210.89000

39、0051488030.320000064986810.930000023756510.850000054592420.480000097047240.920000035801710.600000022157930.140000008500060.400000019573080.880000084805720.500000048692500.130000093819940.160000072500940.230000076880950.520000026321390.190000080209770.290000081746850.020000048766700.320000006824660.4

40、00000096970170.100000071378511.000000066830640.070000064157120.370000009099040.240000096412430.140000071479720.030000097759970.720000067098260.27000003063494 0.70000006738341 0.01000002076566 B =1.0e+002 *2.510000016111824.100000063822142.210000000022821.57000003356329 2.18000002750998 1.87000000445

41、275 * =0.999998307797201.000000225695551.000000193415550.999999093880730.999999968940211.00000066032794x0 =6.181368174725440e-007专业资料整理WORD格式15专业资料整理WORD格式数值分析实验报告x x" x的计算结果为： 6.181368174725440e-007x x（ 2NcondestAcond2AcondA2101.15253088394310232.8905456307542132.89054563075420e+002203.4709596

42、3194066865.5412238417896665.54122384178720e+002306.0505038651128351.1265397557063981.126539755706322e+002e+002e+002403.54948789258247061.3753756968344861.37537569683365e+002506.85501818477940881.1213899375359481.12138993753482e+002601.0820046564093671.7048308151547811.704830815108527e+004e+003e+0037

43、03.2346791451921323.8784811559809363.878481155978439e+003e+002e+002808.31822615391865886.2381429985251386.23814299853018e+002902.0636341434079352.1206963803317052.120696380331079e+003e+002e+0021001.5365928187588971.5591320357384911.559132035738373e+003e+002e+002>> sy5_2please input n:n=8n =8x0

44、 =1.095033343195828e-006x00 =1.705456352162135e-005datx =1.595953017842553e-005给出对x"x x的估计是： 1.705456352162135e-005x xx x" x的理论结果是： 1.095033343195828e-006x x结果相差：1.595953017842553e-005(4)专业资料整理WORD格式16专业资料整理WORD格式数值分析实验报告ncond1Hicond2HicondinfHi4 2.837499999999738 1.551373873892786 2.8374999