数值分析实验报告

1、实验一偏差剖析实验1.1（病态问题）实验目的：算法有“优”与“劣”之分，问题也有“好”与“坏”之别。对数值方法的研究而言，所谓坏问题就是问题自己对扰动敏感者，反之属于好问题。经过本实验可获取一个初步领会。数值剖析的大多数研究课题中，如线性代数方程组、矩阵特点值问题、非线性方程及方程组等都存在病态的问题。病态问题要经过研究和结构特别的算法来解决，自然一般要付出一些代价（如耗用更多的机器时间、占用更多的储存空间等）。问题提出：考虑一个高次的代数多项式明显该多项式的所有根为1,2,20合计20个，且每个根都是单重的。现考虑该多项式的一个扰动此中是一个特别小的数。这相当于是对（1.1）中x19的系数作

2、一个小的扰动。我们希望比较（1.1）和（1.2）根的差异，进而剖析方程（1.1）的解对扰动的敏感性。实验内容：为了实现方便，我们先介绍两个Matlab函数：“roots”和“poly”。此中若变量a储存n+1维的向量，则该函数的输出u为一个n维的向量。设a的元素挨次为a1,a2,an1，则输出u的各重量是多项式方程的所有根；而函数的输出b是一个n+1维变量，它是以n维变量v的各重量为根的多项式的系数。可见“roots”和“poly”是两个互逆的运算函数。上述简单的Matlab程序便获取（1.2）的所有根，程序中的“ess”即是（1.2）中的。实验要求：（1）选择充分小的ess，频频进行上述实验

3、，记录结果的变化并剖析它们。假如扰动项的系数很小，我们自然感觉（1.1）和（1.2）的解应该相差很小。计算中你有什么预料之外的发现表示有些解对于这样的扰动敏感性如何（2）将方程（1.2）中的扰动项改成x18或其余形式，实验中又有如何的现象出现（3）（选作部分）请从理论上剖析产生这一问题的本源。注意我们能够将方程（1.2）写成睁开的形式，同时将方程的解x当作是系数的函数，观察方程的某个解对于的扰动能否敏感，与研究它对于的导数的大小有何关系为何你发现了什么现象，哪些根对于的变化更敏感思虑题一：（上述实验的改良）在上述实验中我们会发现用roots函数求解多项式方程的精度不高，为此你能够考虑用符号函数

4、solve来提升解的精准度，这需要用到将多项式变换为符号多项式的函数poly2sym,函数的详细使用方法可参照Matlab的帮助。实验过程：程序：a=poly(1:20);rr=roots(a);forn=2:21nform=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);y=poly2sym(a);rr=solve(y)forn=2:21nform=1:8ess=10(-6-m);ve=zeros(1,21);ve(n)=ess;

5、a=poly(1:20)+ve;y=poly2sym(a);r=solve(y);-6-ms=max(abs(r-rr)endend数值实验结果及剖析：formatlong-6-m-7-8-9-10n230.9234050060000700008000090000100000110000120000130000140000150000160000170000180000190000200000210000-6-m-11-12-13-14n2030004000050000600007000080000900001000001100001200001300001400001500001600001

6、70000180000190000200000210000议论：利用这种方法进行这种实验，能够很精准的扰动敏感性的一般规律。即当对扰动项的系数愈来愈小时，对其多项式扰动的结果也就愈来愈小，即扰动敏感性与扰动项的系数成正比，扰动项的系数越大，对其根的扰动敏感性就越明显，当扰动的系数一准时，扰动敏感性与扰动的项的幂数成正比，扰动的项的幂数越高，对其根的扰动敏感性就越明显。实验总结：利用MATLAB来进行病态问题的实验，固然其得出的结果是有偏差的，可是可以很简单的得出对一个多次的代数多项式的此中某一项进行很小的扰动，对其多项式的根会有必定的扰动的，因此对于这种病态问题能够借助于MATLAB来进行问题

7、的剖析。学号：06450210姓名：万轩实验二插值法实验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+25x(j)2)i=0,1,n此中l(i)(x),i=0,1,n,n是n次拉格朗日插值基函数。实验要求：选择不

8、停增大的分点数目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)为插值节点结构上述各函数的拉格朗日插值多项式，比较其结果。实验过程：程序：多项式插值的震荡现象（实验2.1）form=1:6subplot(2,3,m)largrang(6*m)%把窗口切割成2*3大小的窗口%对

9、largrang函数进行运转ifm=1title(longn=6)elseifm=2title(longn=12)elseifm=3title(longn=18)elseifm=4title(longn=24)elseifm=5title(longn=30)elseifm=6title(longn=36)end%对每个窗口分别写上标题为插值点的个数end保留为：chazhi.mfunctionlargrang(longn)mm=input(pleaseinputmm(运转第几个函数就输入mm为几):mm=)ifmm=1%d表示定义域的界限值d=1;elseifmm=2|mm=3d=5;endx

10、0=linspace(-d,d,longn);%x的节点ifmm=1y0=1./(1.+25.*x0.2);elseifmm=2y0=x0./(1.+x0.4);elseifmm=3y0=atan(x0);endx=sym(x);n=length(x0);s=0.0;fork=1:np=1.0;forj=1:nifj=kp=p*(x-x0(j)/(x0(k)-x0(j);endends=p*y0(k)+s;endy=s;ifmm=1ezplot(1/(1+25*x2)elseifmm=2ezplot(x/(1+x4)elseifmm=3ezplot(atan(x)endholdonezplot

11、(y,-d,d)holdoff保留为：largrang.m数值实验结果及剖析：对于第一个函数f(x)=1/(1+25x2)对于第二个函数h(x)=x/(1+x4)对于第三个函数g(x)=arctan(x)议论：经过对三个函数得出的largrang插值多项式并在数学软件中的运转，得出函数图象，说了然对函数的支点不是越多越好，而是在函数的两头而言支点越多，而largrang插值多项式不是更为凑近被迫近的函数，反而更为远离函数，在函数两端的跳动性更为明显，argrang插值多项式对函数不收敛。实验总结：利用MATLAB来进行函数的largrang插值多项式问题的实验，固然其得出的结果是有偏差的，可是

12、增添支点的个数进行多次实验，能够找出函数的largrang插值多项式的一般规律，当支点增添时，largrang插值多项式对函数两头不收敛，不是更为迫近，而是更为远离，跳动性更强。因此对于函数的largrang插值多项式问题能够借助于MATLAB来进行问题的剖析，获取比较正确的实验结规律。学号：06450210姓名：万轩实验五解线性方程组的直接方法实验5.1（主元的选用与算法的稳固性）问题提出：Gauss消去法是我们在线性代数中已经熟习的。但因为计算机的数值运算是在一个有限的浮点数会合长进行的，如何才能保证Gauss消去法作为数值算法的稳固性呢Gauss消去法从理论算法到数值算法，其重点是主元的

13、选择。主元的选择从数学理论上看起来平庸，它倒是数值剖析中十分典型的问题。实验内容：考虑线性方程组编制一个能自动选用主元，又好手动选用主元的求解线性方程组的Gauss消去过程。实验要求：61786115（1）取矩阵A,b，则方程有解x*(1,1,1)T。取n=10861158614计算矩阵的条件数。让程序自动选用主元，结果如何2）现选择程序中手动选用主元的功能。每步消去过程总选用按模最小或按模尽可能小的元素作为主元，察看并记录计算结果。若每步消去过程总选用按模最大的元素作为主元，结果又如何剖析实验的结果。3）取矩阵阶数n=20或许更大，重复上述实验过程，察看记录并剖析不一样的问题及消去过程中选择

14、不一样的主元时计算结果的差异，说明主元素的选用在消去过程中的作用。4）选用其余你感兴趣的问题或许随机生成矩阵，计算其条件数。重复上述实验，察看记录并剖析实验结果。实验过程：程序：成立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)pausem,n=size(A);nb=n+1;Ab=Abr=input(请输入能否为手动，手动输入1，自

15、动输入0：r=)fori=1:n-1ifr=0pivot,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(

16、n,n);fori=n-1:-1:1x(i)=(Ab(i,nb)-Ab(i,i+1:n)*x(i+1:n)/Ab(i,i);end数值实验结果及剖析：取矩阵A的阶数:n=10，自动选用主元：formatlonggauss请输入矩阵A的阶数:n=10n=10条件数对应的范数是p-范数：p=1p=1请输入能否为手动，手动输入1，自动输入0：r=0r=0取矩阵A的阶数:n=10，手动选用主元：选用绝对值最大的元素为主元：gauss请输入矩阵A的阶数:n=10n=10条件数对应的范数是p-范数：p=2p=2pp=请输入能否为手动，手动输入1，自动输入0：r=1r=1ans=1111111111选用绝对

17、值最小的元素为主元：gauss请输入矩阵A的阶数:n=10n=10条件数对应的范数是p-范数：p=2p=2pp=03请输入能否为手动，手动输入1，自动输入0：r=1r=1ans=1.000000000000001.000000000000001.000000000000011.00000000000003取矩阵A的阶数:n=20，手动选用主元：选用绝对值最大的元素为主元：gauss请输入矩阵A的阶数:n=20条件数对应的范数是p-范数：p=1p=1pp=ans=选用绝对值最小的元素为主元：gauss请输入矩阵A的阶数:n=20.n=20条件数对应的范数是p-范数：p=2p=2pp=请输入能否为

18、手动，手动输入1，自动输入0：r=1r=1ans=1.000000000000001.000000000000001.000000000000011.00000000000006999891.000000000000231.000000000000901.000000000003521.000000000012731.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

19、=请输入能否为手动，手动输入1，自动输入0：r=1r=1ans=1.000000000000511.000000000313541.000000002688051.00000000084337gauss请输入矩阵A的阶数:n=7n=7条件数对应的范数是p-范数：p=2p=2pp=请输入能否为手动，手动输入1，自动输入0：r=1r=1ans=1.000000000043371.000000001211431.00000000152825该问题在主元选用与算出结果有着很大的关系，取绝对值大的元素作为主元比取绝对值小的元素作为主元时产生的结果比较正确，即选用绝对值小的主元时结果产生了较大的偏差，条件

20、数越大产生的偏差就越大。议论：在gauss消去法解线性方程组时，主元的选择与算法的稳固性有亲密的联系，选用绝对值大的元素作为主元比绝对值小的元素作为主元时对结果产生的偏差较小。条件数越大对用gauss消去法解线性方程组时，对结果产生的偏差就越大。实验总结：对用gauss消去法解线性方程组时，主元的选用与算法的稳固性有亲密的联系，选用合适的主元有益于得出稳固的算法，在算法的过程中，选用绝对值较大的主元比选用绝对值较小的主元更有益于算法的稳固，选用绝对值最大的元素作为主元时，得出的结果相对较正确较稳固。条件数越小，对用这种方法得出的结果更正确。在算除法的过程中要尽量防止使用较小的数做为除数，免得发

21、生结果数目级加大，使大数吃掉小数，产生舍入偏差。学号：06450210姓名：万轩实验5.2（线性代数方程组的性态与条件数的预计）问题提出：理论上，线性代数方程组Ax矩阵的条件数的确是对矩阵病态性的刻画，b的摄动知足但在实质应用中直接计算它明显不现实，因为计算A1往常要比求解方程Axb还困难。实验内容：Matlab中供给有函数“condest”能够用来预计矩阵的条件数，它给出的是按1-范数的条件数。第一结构非奇怪矩阵A和右端，使得方程是能够精确求解的。再人为地引进系数矩阵和右端的摄动A和b，使得A和b充分小。实验要求：（1）假定方程Ax=b的解为x，求解方程(AA)x?bb，以1-范数，给出x?

22、x的计算结果。x2）选择一系列维数递加的矩阵（能够是随机生成的），比较函数“condest”所需机器时间的差异.考虑若干逆是已知的矩阵，借助函数“eig”很简单给出cond2(A)的数值。将它与函数“cond(A,2)”所获取的结果进行比较。（3）利用“condest”给出矩阵A条件数的预计，针对（1）中的结果给出x的x理论预计，并将它与（1）给出的计算结果进行比较，剖析所得结果。注意，假如给出了cond(A)和A的预计，立刻就能够给出A1的预计。（4）预计着名的Hilbert矩阵的条件数。实验过程：程序：n=input(pleaseinputn:n=)%输入矩阵的阶数a=fix(100*ra

23、nd(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)%得出?x的理论结果xxxx保留为：fanshu.mfunctionx=geshow(A,B)%用高斯消去法解方程组m,n=size(A);nb=n+1;AB=AB;fori=1:n-1pivot=AB(i,i);for

24、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);fori=n-1:-1:1x(i)=(AB(i,nb)-AB(i,i+1:n)*x(i+1:n)/AB(i,i);end保留为：geshow.mfunctioncond2(A)B=A*A;V1,D1=eig(B);V2,D2=eig(B(-1);%自定义求二阶条件数cond2A=sqrt(max(max(D1)*sqrt(max(max(D2)end保留为：cond2.mformatlongforn=1

25、0:10:100n=n%n为矩阵的阶A=fix(100*randn(n);%随机生成矩阵AcondestA=condest(A)%用condest求条件数cond2(A)%用自定义的求条件数condA2=cond(A,2)%用cond求条件数pause%运转一次暂停end保留为：shiyan52.mn=input(pleaseinputn:n=)%输入矩阵的阶数a=fix(100*rand(n)+1;%随机生成一个矩阵ax=ones(n,1);%假定知道方程组的解全为1b=a*x;%用矩阵a和以知解得出矩阵bdata=rand(n)*0.00001;%随即生成扰动矩阵datadatb=rand

26、(n,1)*0.00001;%随即生成扰动矩阵datbA=a+data;B=b+datb;xx=geshow(A,B);%利用第一小问的geshow.m求出解阵x0=norm(xx-x,1)/norm(x,1)xx?x%得出的理论结果xxx00=cond(A)/(1-norm(inv(A)*norm(xx-x)*(norm(xx-x)/(norm(A)+norm(datb)/norm(B)%得出x?xx的预计值xxdatx=abs(x0-x00)%求二者之间的偏差保留为：sy5_2.mformatlongforn=4:11n=n%n为矩阵的阶数Hi=hilb(n);%生成Hilbert矩阵co

27、nd1Hi=cond(Hi,1)%求Hilbert矩阵得三种条件数cond2Hi=cond(Hi,2)condinfHi=cond(Hi,inf)pauseend数值实验结果及剖析：fanshupleaseinputn:n=6n=6a=142516881989329385489260144088501316235219292324010100737241437227701x=111111b=251410221157218187data=1.0e-005*datb=1.0e-005*A=1.0e+002*B=1.0e+002*=x0=x?x的计算结果为：x2）NcondestAcond2Acon

28、dA210e+00220e+00230e+002e+002e+00240e+00250e+00260e+004e+003e+00370e+003e+002e+00280e+00290e+003e+002e+002100e+003e+002e+002sy5_2pleaseinputn:n=8n=8给出对x?的预计是：xxxx?x的理论结果是：x结果相差：(4)ncond1Hicond2HicondinfHi4e+004e+004e+0045e+005e+005e+0056e+007e+007e+0077e+008e+008e+0088e+010e+010e+0109e+012e+011e+012

29、10e+013e+013e+01311e+015e+014e+015议论：线性代数方程组的性态与条件数有着很重要的关系，既矩阵的条件数是刻画矩阵性质的一个重要的依照，条件数越大，矩阵“病态”性越严重，在解线性代数方程组的过程中较简单产生比较大的偏差，则在实质问题的操作过程中，我们一定要减少对条件数来求解，把条件数较大的矩阵化成条件数较小的矩阵来进行求解。实验总结：在本次实验中，使我们知道了矩阵条件数对线性代数方程组求解的影响，条件数越大，对最后解的影响的越大，hilbert矩阵是一个很”病态”的矩阵，他的条件数跟着阶数的增添而增大，每增添一阶，条件数就增大一个数目级，在求解的过程中要尽量防止h

30、ilbert矩阵学号：06450210姓名：万轩实验七非线性方程求根实验7.1（迭代法、初始值与收敛性）实验目的：初步认识非线性问题的迭代法与线性问题迭代法的差异，商讨迭代法及初始值与迭代收敛性的关系。问题提出：迭代法是求解非线性方程的基本思想方法，与线性方程的状况同样，其结构方法能够有多种多样，但重点是如何才能使迭代收敛且有较快的收敛速度。实验内容：考虑一个简单的代数方程针对上述方程，能够结构多种迭代法，如在实轴上取初始值x0，请分别用迭代（7.1）-（7.3）作实验，记录各算法的迭代过程。实验要求：1）取定某个初始值，分别计算（7.1）-（7.3）迭代结果，它们的收敛性如何重复选用不一样的

31、初始值，频频实验。请自选设计一种比较形象的记录方式（如利用Matlab的图形功能），剖析三种迭代法的收敛性与初值选用的关系。2）对三个迭代法中的某个，取不一样的初始值进行迭代，结果如何试剖析迭代法对不一样的初值能否有差异3）线性方程组迭代法的收敛性是不依靠初始值选用的。比较线性与非线性问题迭代的差异，有何结论和问题。实验过程：程序：clearclcs=input(请输入要运转的方程，运转第几个输入几s=);clfifs=1%决定坐标轴的范围和初始值a=-1.5;b=2.5;y00=0;x00=input(请输入第一个函数的初值：x00=);elseifs=2a=0.1;b=6.5;y00=0;

32、x00=input(请输入第二个函数的初值：x00=);elseifs=3a=0;b=2;y00=0;x00=input(endx=linspace(a,b,80);y0=x;y1=zxy7f(x,s);cleary;y=y0;y1;ifs=1plot(x,y,linewidth,1)legend(y=x,y=f1)title(x(n+1)=x(n)2-1)elseifs=2plot(x,y,linewidth,2)legend(y=x,y=f2)请输入第三个函数的初值：x00=);%计算直线y=x计算迭代函数y=f(x)%绘图输出标题title(x(n+1)=1+1/x(n)elseifs=

33、3plot(x,y,linewidth,3)legend(y=x,y=f3)title(x(n+1)=sqrtx(n)+1)endholdonplot(ab,0,0,k-,00,ab,k-)axis(a,b,a,b)%画坐标轴z=;fori=1:15%画蛛网图，迭代过程为n=15次xt(1)=x00;yt(1)=y00;%决定始点坐标xt(2)=zxy7f(xt(1),s);%决定终点坐标yt(2)=zxy7f(xt(1),s);zxyplot7(xt,yt,0.6)%画蛛网图ifi=5pause%按随意键逐次察看前5次迭代的蛛网图endx00=xt(2);y00=yt(2);z=z,xt(1);%将本次迭代的终点作为下次的始点%保留迭代点end保留为：zxy7.mfunctiony=zxy7f(x,s)ifs=1y=(x.*x-1);elseifs=2y=(1+1./x);elseifs=3y=sqrt(x+1);end保留为：zxy7f.mfunctionout=zxyplot7(x,y,p)u(1)=0;v(1)=(y(2)-y(1);%画一次迭代的蛛网图，改变p调理箭头的大小%画出始点(x(1),y(1)终点(x(2),y(2)的有向折线段u(2)=eps;v(2)=eps;h=quiver(x(1)x(1),y(1)y(2),u,v