数值分析实验报告

数值分析实验报告姓名：学号：学院：实验一病态问题实验目的：算法有“优”与“劣”之分，问题也有“好”与“坏”之别。对数值方法的研究而言，所谓坏问题就是问题本身对扰动敏感者，反之属于好问题。通过本实验可获得一个初步体会。数值分析的大部分研究课题中，如线性代数方程组、矩阵特征值问题、非线性方程及方程组等都存在病态的问题。病态问题要通过研究和构造特殊的算法来解决，当然一般要付出一些代价（如耗用更多的机器时间、占用更多的存储空间等）。问题提出：考虑一个高次的代数多项式201()(1)(2)(20)()kpxxxxxk（E.1.1）显然该多项式的全部根为1,2,…,20共计20个，且每个根都是单重的。现考虑该多项式的一个扰动19()0pxx(E.1.2)其中是一个非常小的数。这相当于是对（E.1.1）中19x的系数作一个小的扰动。我们希望比较（E.1.1）和（E.1.2）根的差别，从而分析方程（E.1.1）的解对扰动的敏感性。实验内容：为了实现方便，我们先介绍两个Matlab函数：“roots”和“poly”。roots(a)u其中若变量a存储n+1维的向量，则该函数的输出u为一个n维的向量。设a的元素依次为121,,,naaa，则输出u的各分量是多项式方程01121nnnnaxaxaxa的全部根；而函数bpolyv的输出b是一个n+1维变量，它是以n维变量v的各分量为根的多项式的系数。可见“roots”和“poly”是两个互逆的运算函数。);21,1(zerosve;)2(essve))20:1((vepolyroots上述简单的Matlab程序便得到（E.1.2）的全部根，程序中的“ess”即是（E.1.2）中的。实验要求：（1）选择充分小的ess，反复进行上述实验，记录结果的变化并分析它们。如果扰动项的系数很小，我们自然感觉（E.1.1）和（E.1.2）的解应当相差很小。计算中你有什么出乎意料的发现？表明有些解关于如此的扰动敏感性如何？（2）将方程（1.2）中的扰动项改成18x或其它形式，实验中又有怎样的现象出现？实验过程：根据上述方法，为了分析所有解关于扰动的敏感性，实验中在一定范围内改变扰动项系数以及改变扰动项，对扰动前后的解的值进行对比，观察扰动结果。实验程序如下：程序：a=poly(1:20);r1=roots(a);forn=2:21form=1:9ess=10^(-6-m);%改变扰动系数的大小ve=zeros(1,21);ve(n)=ess;%改变扰动项的位置r2=roots(a+ve);s=max(abs(r2-r1))%计算变化前后变化最大的解的差值endend数值实验结果及分析：表1.1改变扰动项系数值和扰动项前后解的变化情况-6-mn-7-8-9-1022.797226874783311.867536320091581.060527623807480.2527314421904731.693766997674240.923106667069640.084716145697410.4080402640941140.854013934155360.199410220200610.03972935295834050.110311005388710.042965323628440060000700008000090000100000110000120000130000140000150000160000170000180000190000200000210000-6-mn-11-12-13-1420.038776764393800.162565848682800.13322664013598030.02164258317546000400005000060000700008000090000100000110000120000130000140000150000160000170000180000190000200000210000结果分析：根据表1.1所示的实验结果，我们可以很精确的看到扰动敏感性的一般规律。即当对扰动项的系数越来越小时，对其多项式扰动的结果也就越来越小，即扰动敏感性与扰动项的系数成正比，扰动项的系数越大，对其根的扰动敏感性就越明显，当扰动的系数一定时，扰动敏感性与扰动的项的幂数成正比，扰动的项的幂数越高，对其根的扰动敏感性就越明显。实验总结：利用MATLAB来进行病态问题的实验，虽然其得出的结果是有误差的，但是我们还是可以看出一些规律，对一个多次的代数多项式的其中某一项进行很小的扰动，对其多项式的根会有一定的扰动的，而且扰动项的幂数越高，则影响越大，借助于MATLAB我们很容易对这类病态问题进行问题的分析。实验二多项式插值的振荡现象问题提出：考虑一个固定的区间上用插值逼近一个函数。显然拉格朗日插值中使用的节点越多，插值多项式的次数就越高。我们自然关心插值多项式的次数增加时，nLx是否也更加靠近被逼近的函数。Runge给出了一个极著名例子。设区间[-1，1]上函数21.125fxx实验内容：考虑区间[-1，1]的一个等距划分，分点为：21,0,1,2,,,iixinn则拉格朗日插值多项式为：201.125nniiiLxlxx其中，,0,1,2,,ilxin是n次Lagrange插值基函数。实验要求：⑴选择不断增大的分点数目n=2,3…，画出原函数fx及插值多项式函数nLx在[-1，1]上的图像，比较分析实验结果。（2）选择其它的函数，例如定义在区间[-5，5]上的函数4,arctan,1xhxgxxx重复上述的实验看其结果如何。实验过程：程序：form=1:6subplot(2,3,m)%把窗口分割成2*3大小的窗口largrang(4*m)%对largrang函数进行运行ifm==1title('n=4')elseifm==2title('n=8')elseifm==3title('n=12')elseifm==4title('n=16')elseifm==5title('n=20')elseifm==6title('n=24')end%对每个窗口分别写上标题为插值点的个数endfunctionlargrang(longn)mm=input('pleaseinputmm(运行第几个函数就输入mm为几):mm=')ifmm==1%d表示定义域的边界值d=1;elseifmm==2||mm==3d=5;endx0=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*x^2)')elseifmm==2ezplot('x/(1+x^4)')elseifmm==3ezplot('atan(x)')endholdonezplot(y,[-d,d])holdoff数值实验结果及分析：图2.1函数21125fxx的多项式差值图2.2函数41xhxx的多项式差值-1-0.500.510.050.10.150.20.250.3xn=4-1-0.500.5100.20.40.60.8xn=8-1-0.500.5100.20.40.60.81xn=12-1-0.500.5100.511.5xn=16-1-0.500.51-0.500.511.5xn=20-1-0.500.51-0.500.511.5xn=24-505-0.6-0.4-0.200.20.40.6xn=6-505-1.5-1-0.500.511.5xn=12-505-1.5-1-0.500.511.5xn=18-505-4-2024xn=24-505-3-2-1012xn=30-505-30-20-1001020xn=36图2.3函数arctangxx的多项式差值由上述三个图，对于第一个函数，试验中插值的点数为4,8,…,24,对第二，三个函数，插值点数均为6,12，…,36,用插值法和函数本身图像进行对比图如上述所示。通过对三个函数得出的largrang插值多项式得出函数图像，说明了对函数的支点不是越多越好，而是在函数的两端而言支点越多，而largrang插值多项式不是更加靠近被逼近的函数，反而更加远离函数，在函数两端的跳动性更加明显，largrang插值多项式对函数不收敛。实验总结：利用MATLAB来进行函数的largrang插值多项式问题的实验，虽然其得出的结果是有误差的，但是增加支点的个数进行多次实验，可以找出函数的largrang插值多项式的一般规律，当支点增加时，largrang插值多项式对函数两端不收敛，不是更加逼近，而是更加远离，跳动性更强。所以对于函数的largrang插值多项式问题可以借助于MATLAB来进行问题的分析，得到比较准确的实验结规律。-505-1.5-1-0.500.511.5xn=6-505-2-1012xn=12-505-3-2-10123xn=18-505-3-2-10123xn=24-505-3-2-10123xn=30-505-3-2-10123xn=36实验三线性代数方程组的性态与条件数的估计问题提出：理论上，线性代数方程组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,)(,,实验过程：程序：⑴n=input('pleaseinputn:n=')；%输入矩阵的阶数a=fix(100*rand(n))+1；%随机生成一个矩阵ax=ones(n,1)；b=a*x；%用矩阵a和以知解得出矩阵bdata=rand(n)*0.00001；%随即生成扰动矩阵datadatb=rand(n,1)*0.00001；%随即生成扰动矩阵datbA=a+data；B=b+datb；xx=geshow(A,B)；x0=norm(xx-x,1)/norm(x,1)；%得出xxxxxˆ的理论结果%%用高斯消去法解方程组%%functionx=geshow(A,B)[m,n]=size(A);nb=n+1;AB=[AB];fori=1:n-1pivot=AB(i,i);fork=i+1:nAB(k,i:nb)=AB(k,i:nb)-(AB(k,i)/pivot)*AB(i,i:nb);endendx=zeros