数值分析实验报告

《数值分析》实验报告1实验1.1病态问题1.1程序代码%对应于实验1.1的程序clear;prompt={'请输入扰动项n:0~20间的整数','请输入扰动系数ess:0～1间的小数'};def={'19','0.00001'};dlgTitle='实验1.1病态问题';lineNo=1;result=inputdlg(prompt,dlgTitle,lineNo,def);NE=str2num(char(result));%转换结果if((NE(1,1)20)|(NE(1,1)0)),error('请输入合适的项数！');endess=NE(2,1);VE=zeros(1,21);VE(21-NE(1,1))=ess;%将扰动项变成向量模式，以便相加ROOT=roots(poly(1:20)+VE);%求p(x)+ess*x=0;的解1.2结果记录1.2.1在第19项加扰动ess后，方程0)(19201xesskxk(1)，的根如表1所示：阶数\ess1.00E-011.00E-031.00E-051.00E-071.00E-09136.4232+13.18i27.0817+5.03812i22.5961+2.3083i20.422+0.999203i19.952+0i236.4232-13.18i27.0817-5.03812i22.5961-2.3083i20.422-0.999203i19.2293+0i318.2486+16.1526i19.5337+9.1664i18.8972+5.00563i18.1572+2.4702i17.6573+0.692896i418.2486-16.1526i19.5337-9.1664i18.8972-5.00563i18.1572-2.4702i17.6573-0.692896i511.3192+10.81i13.8235+7.77167i14.9123+4.95848i15.3149+2.69865i15.4524+0.875524i611.3192-10.81i13.8235-7.77167i14.9123-4.95848i15.3149-2.69865i15.4524-0.875524i78.773+6.95838i10.7211+5.4609i12.0289+3.73551i12.8466+2.06246i13.3527+0.486992i88.773-6.95838i10.7211-5.4609i12.0289-3.73551i12.8466-2.06246i13.3527-0.486992i97.47593+4.38779i8.91282+3.47317i10.059+2.33021i10.9216+1.10366i11.8578+0i107.47593-4.38779i8.91282-3.47317i10.059-2.33021i10.9216-1.10366i11.0427+0i116.60562+2.56677i7.69268+1.89884i8.63828+1.0564i9.56629+0i9.9916+0i126.60562-2.56677i7.69268-1.89884i8.63828-1.0564i9.11508+0i9.00201+0i135.90709+1.20368i6.75761+0.654714i7.70896+0i7.99387+0i7.99952+0i145.90709-1.20368i6.75761-0.654714i7.028+0i7.00027+0i7.00009+0i155.291225.95208+0i5.99942+0i66165.103695.00061+0i5.00001+0i55173.999784444183333319222222011111表11.2.2在第18项加扰动ess后，方程0)(18201xesskxk(2)，的根如表2所示：阶数\ess1.00E-011.00E-031.00E-051.00E-071.00E-09129.1493+5.98286i23.724+2.8203i20.994+1.34421i19.7656+0.300013i19.9983229.1493-5.98286i23.724-2.8203i20.994-1.34421i19.7656-0.300013i19.0115320.2125+11.071i19.4236+6.14475i18.4909+3.24042i17.8435+1.30236i17.9649420.2125-11.071i19.4236-6.14475i18.4909-3.24042i17.8435-1.30236i17.0524513.4752+9.28466i14.8344+6.09209i15.3107+3.55896i15.4553+1.5856i15.9662613.4752-9.28466i14.8344-6.09209i15.3107-3.55896i15.4553-1.5856i14.9771710.1007+6.44431i11.6561+4.62319i12.6114+2.82705i13.1902+1.20025i14.0752810.1007-6.44431i11.6561-4.62319i12.6114-2.82705i13.1902-1.20025i12.904798.29727+4.13106i9.60113+2.99431i10.579+1.73097i11.275+0.452354i12.0762108.29727-4.13106i9.60113-2.99431i10.579-1.73097i11.275-0.452354i10.961117.14293+2.37702i8.18209+1.56759i9.04344+0.625845i9.93183+0i10.0166127.14293-2.37702i8.18209-1.56759i9.04344-0.625845i9.00952+0i8.99504136.27379+1.02957i7.0834+0.388249i7.93758+0i7.99935+0i8.00115146.27379-1.02957i7.0834-0.388249i7.00355+0i7+0i6.9998155.68402+0i6.000026.000026.000026.00002165.0129455551744444183333319222222011111表21.3结果分析多项式201)(kkx展开后，从高到低，其各次项的系数依次如下：1-21020615-1.2569e+0065.3328e+007-1.6723e94.0172e10-7.5611e111.131e13-1.3559e+0141.3075e15-1.0142e166.3031e16-3.1133e+0171.2066e+018-3.6e+0188.0378e+018-1.2871e+0191.3804e+019-8.7529e182.4329e181.3.1对于方程0)(201kkx（3），其方程左边的多项式的第19次项的系数为-210，在其上加上ess（1e-9~0.1）的扰动，扰动的程度不大于0.05％，但是却给方程(3)的解带来了很大的颠覆性变化！详细分析发现，在第19次项的系数上的加ess后，对第20、19、18、17次的对应解无影响；但是对其后各次的解，影响却是愈来愈大，直至完全颠覆。如x1原来的解为20，加入ess=0.1后，变为29.1493+5.98286i。1.3.2对于方程(3)，其左边多项式的第18次项的系数为20615，在其上加上ess（1e-9~0.1）的扰动。分析发现，对第20、19、18、17、16次的对应解无影响；但是对其后各次的解，影响却是颠覆性的。这与在第19项加ess的结果是一致的。只是程度没有前者大。1.3.3问题的根源求解实系数高次多项式方程，并没有通用的代数解法，而只能通过迭代法求解。这就为上述的不稳定现象奠定了伏笔。同时，对方程（1）或（2）两边对ess求导（认定x是关于ess的函数）发现，由于越低次项的系数越大，故越低次项对应的解对dx/dess的波动就越敏感！这也就解释了为什么无论是扰动项在18还是19，前5个高次项对应的解基本不变，而后面的却大受其害的现象。2实验2.1多项式插值的振荡现象2.1程序代码%对应于实验2.1的程序clear;prompt={'请输入实验函数:f(x)或和h(x)的代号，前者为f，后者是h'};def={'f'};dlgTitle='实验2.1多项式插值的振荡现象';lineNo=1;result=inputdlg(prompt,dlgTitle,lineNo,def);fun=char(result);%转换结果if((fun~='f')&(fun~='h')),error('请选择合适的函数！');break;endif(fun=='f')a=-1;b=1;f=inline('1./(1+25*x.^2)');elsea=-5;b=5;f=inline('x./(1+4*x.^4)');endfplot(f,[ab],'m');fori=2:2:10holdon;x0=linspace(a,b,i+1);y0=feval(f,x0);x=a:0.1:b;y=Lagrange(x0,y0,x);if(i==2)plot(x,y,'gp');plot(x,y,'g:');elseif(i==4)plot(x,y,'r+');plot(x,y,'r-.');elseif(i==6)plot(x,y,'kd');plot(x,y,'k-');elseif(i==8)plot(x,y,'co');plot(x,y,'c--');elseplot(x,y,'b*');plot(x,y,'b:');end;end;2.2结果记录函数f(x)=1/(1+25*x^2)及其2、4、6、8、10点插值函数在[-1,1]上的图象：图2－1函数h(x)=x/(1+4*x^4)及其2、4、6、8、10点插值函数在[-1,1]上的图象：图2－22.3结果分析从上述图2－1和图2－2不难发现插值点数从2点增到4点，插值函数对原函数的近似程度要加大。但是，继续增多插值点后，如6、8、10点，插值函数对原函数的逼近效果并没有相应提高，反而出现了振荡现象！3实验3.1编制多项式最小二乘拟合函数3.1程序代码clear;prompt={'请输入进行最小二乘拟合的多项式的次数n：','请输入Xi','请输入Yi'};def={'3','-1:0.5:2','[-4.447-0.4520.5510.048-0.4470.5494.522]'};dlgTitle='实验3.1编制插值函数，并进行曲线拟合';lineNo=1;result=inputdlg(prompt,dlgTitle,lineNo,def);n=str2num(char(result(1,1)));%转换结果x0=str2num(char(result(2,1)));y0=str2num(char(result(3,1)));if((n1)|(n9)),error('请选择合适的次数！');break;end;if(length(x0)~=length(y0))error('Xi序列和Yi序列的个数必须一致！');break;end;plot(x0,y0,'rp');%输出原数据gridon;alph=polyfit(x0,y0,n);%利用库函数计算进行最小二乘拟合的n次多项式的各项系数。x=-1:0.01:2;y=polyval(alph,x);%利用库函数求取拟合函数值holdon;plot(x,y);y=polyval(alph,x0);%利用库函数求取拟