《数值逼近课程设计.doc》由会员分享，可在线阅读，更多相关《数值逼近课程设计.doc（11页珍藏版）》请在三一办公上搜索。
1、第一题 龙格现象1901年Runge给出一个例子 定义在区间-1,1上。这是一个很光滑的函数，它的任意阶导数都存在。对它在-1,1上做等距节点插值时，会发现插值多项式在变量靠近-1和1时余项会随着n的增大而很大；另一个现象就是插值多项式随节点增多而振动更多。这种插值多项式当节点增多时反而不能更好的接近被插函数的现象称为Runge现象。function f=runge(n)syms t z z1; x=-1:2/n:1;y=1./(1+25*x.2);z=0;for i=1:(n+1) z1=y(1,i); for j=1:(n+1) if(i=j) z1=z1.*(t-x(1,j)./(x(1
2、,i)-x(1,j); end end z=z+z1;endf=expand(z);在命令窗口输入命令：y1=runge(4)y1=1250/377*t4-3225/754*t2+1; y2=runge(8)y2=1+228601250/3725137*t4-383000000/3725137*t6+200000000/3725137*t8-98366225/7450274*t2; y3=runge(12)y3=1+367051586875/1847048164*t4-107641853578125/112669938004*t6+62017871484375/28167484501*t8+2
3、5628906250000/28167484501*t12-551599221900/28167484501*t2-65809335937500/28167484501*t10;x=-1:0.01:1;y=(1+25*x.2).(-1);t=-1:0.01:1;y1=1250/377.*t.4-3225/754.*t.2+1;y2=1+228601250/3725137.*t.4-383000000/3725137.*t.6+200000000/3725137.*t.8-98366225/7450274.*t.2;y3=1+367051586875/1847048164.*t.4-107641
4、853578125/112669938004.*t.6+62017871484375/28167484501.*t.8+25628906250000/28167484501.*t.12-551599221900/28167484501.*t.2-65809335937500/28167484501.*t.10;plot(x,y,t,y1,:,t,y2,-.,t,y3,*)legend(f(x),n=4,n=8,n=12)图1结论： 由上图可以看出用高次插值多项式是不妥当的，误差会随着阶数的提高越来越严重。因此，为提高插值的精度，可以采用分段插值。第二题 Chebyshev多项式画出Chebys
5、hev多项式在n=2，3，4，5时的图像x=-1:0.1:1;T2=cos(2*acos(x);T3=cos(3*acos(x);T4=cos(4*acos(x);T5=cos(5*acos(x);plot(x,T2,x,T3,-.,x,T4,:,x,T5,-)legend(n=2,n=3,n=4,n=5)图2第三题 Remez算法1.求函数f(x)=在-1,1上的二次多项式逼近。编写M文件function remezs(F,x1,x2,x3,x4,e,min,max)format long;syms x ;y1=F;xA=5;xB=5;xC=5;xD=5;m=0;n=0;X1(1)=x1;X
6、2(1)=x2;X3(1)=x3;X4(1)=x4;P1(1)=subs(F,x,x1);P2(1)=subs(F,x,x2);P3(1)=subs(F,x,x3);P4(1)=subs(F,x,x4); d=2; while(abs(abs(subs(y1,x,xA)-abs(m)e | abs(abs(subs(y1,x,xB)-abs(m)e | abs(abs(subs(y1,x,xC)-abs(m)e | abs(abs(subs(y1,x,xD)-abs(m)e) a=1,1,x1,x12;1,-1,x2,x22;1,1,x3,x32;1,-1,x4,x42; z1=subs(F,
7、x,x1); z2=subs(F,x,x2); z3=subs(F,x,x3); z4=subs(F,x,x4); b=z1;z2;z3;z4; g=ab; c=g(4,1);g(3,1);g(1,1); m=g(2,1); p=poly2sym(c); P1(d)=subs(p,x,x1);P2(d)=subs(p,x,x2);P3(d)=subs(p,x,x3);P4(d)=subs(p,x,x4); y1=F-p; y2=diff(y1); y3=diff(y2); s=(max-min)/50; j=1; for i=s:s:max-min k1=subs(y2,x,min+i-s);
8、k2=subs(y2,x,min+i); if(k10)|(k10)&(k20.00000001 | abs(abs(xC)-abs(xc)0.00000001) y4=xb-y2/y3; xB=subs(y4,x,xb); t=xb; xb=xB; xB=t; y5=xc-y2/y3; xC=subs(y5,x,xc); t=xc; xc=xC; xC=t;end x1=xA;X1(d)=x1; x2=xB;X2(d)=x2; x3=xC;X3(d)=x3; x4=xD;X4(d)=x4; n=n+1;d=d+1; end nmcp=poly2sym(c);y=p X1 X2 X3 X4 P
9、1 P2 P3 P4 在命令窗口输入命令：syms xy=exp(x);remezs (y,-1,-0.5,0.5,1,0.000000001,-1,1) %当初始点组为-1,-0.5,0.5,1时n = 3 %迭代步数m = -0.04501738840272 %最佳逼近c = 0.55404090635670 %最佳逼近多项式系数 1.13018380524108 0.98903972845855y=1247589209708019/2251799813685248*x2+5089895364143919/4503599627370496*x+8908477905081045/900719
10、9254740992 %最佳逼近多项式X1 = -1 -1 -1 -1X2 = -0.500000000000000 -0.438621271668792 -0.436958142528573 -0.436958064363172X3 = 0.500000000000000 0.560938664934263 0.560058700967099 0.560057761761697X4 = 1 1 1 1P1 = 0.367879441171442 0.412216302056878 0.412896544738883 0.412896829574159P2 = 0.60653065971263
11、3 0.562193798827198 0.599907881178971 0.600981082304536P3 = 1.648721270700128 1.693058131585564 1.797333670243964 1.795792657883955P4 = 2.718281828459046 2.673944967573610 2.673264724891605 2.673264440056329remezs (y,-1,0,0.5,1,0.000000001,-1,1) %当初始点组为-1,0,0.5,1时n = 4 %迭代步数m = -0.04501738840244 %最佳
12、逼近c = 0.55404090635629 %最佳逼近多项式系数 1.13018380524136 0.98903972845896y=4990356838828379/9007199254740992*x2+2544947682072583/2251799813685248*x+8908477905084741/9007199254740992 %最佳逼近多项式X1 = -1 -1 -1 -1 -1X2 = Columns 1 through 4 0 -0.421854446335344 -0.439110722460123 -0.436957524811125 Column 5 -0.4
13、36958064364547X3 = Columns 1 through 4 0.500000000000000 0.616666867105925 0.559560192587622 0.560059523299712 Column 5 0.560057761761412X4 = 1 1 1 1 1P1 = Columns 1 through 4 0.367879441171442 0.401056989982813 0.412479164166748 0.412896473452533 Column 5 0.412896829573882P2 = Columns 1 through 4 1
14、.000000000000000 0.966822451188630 0.611229767837094 0.599592370658370 Column 5 0.600981481349462P3 = Columns 1 through 4 1.648721270700128 1.681898819511499 1.897342025310450 1.794919743126876 Column 5 1.795794097603871P4 = Columns 1 through 4 2.718281828459046 2.685104279647675 2.673682105463739 2
15、.673264796177954 Column 5 2.673264440056606remezs (y,-0.8,-0.3,0.4,0.7,0.000000001,-1,1) %当初始点组为-0.8,-0.3,0.4,0.7时n = 5 %迭代步数m =-0.04501738840282 %最佳逼近c = 0.55404090635688 %最佳逼近多项式系数 1.13018380524098 0.98903972845837y=2495178419416851/4503599627370496*x2+5089895364143459/4503599627370496*x+890847790
16、5079421/9007199254740992 %最佳逼近多项式X1 = -0.800000000000000 -1.000000000000000 -1.000000000000000 -1.000000000000000 -1.000000000000000 -1.000000000000000X2 = -0.300000000000000 -0.389100143282599 -0.456582740283911 -0.436909973732124 -0.436958078197744 -0.436958064362609X3 = 0.400000000000000 0.363012
17、301238986 0.556102082418996 0.560206774246383 0.560057758708891 0.560057761761852X4 = 0.700000000000000 1.000000000000000 1.000000000000000 1.000000000000000 1.000000000000000 1.000000000000000P1 = 0.449328964117222 0.466303115476679 0.408989321350584 0.412867342011747 0.412896826843031 0.4128968295
18、74262P2 = 0.740818220681718 0.723844069322260 0.636556522698576 0.588456694999194 0.601012202753950 0.600981123862051P3 = 1.491824697641270 1.508798849000728 1.478763424231403 1.788849706405623 1.796051917382699 1.795791008202466P4 = 2.013752707470477 1.996778556111019 2.677171948279904 2.6732939276
19、18741 2.673264442787456 2.6732644400562272.画出f(x)=函数在中进行二次多项式最佳逼近逼近，Remez算法中交错点组随迭代次数k的变化规律，以及逼近误差曲线的形象syms xx1=-1;x2=-0.5;x3=0.5;x4=1;F=exp(x);a=1,1,x1,x12;1,-1,x2,x22;1,1,x3,x32;1,-1,x4,x42;z1=subs(F,x,x1);z2=subs(F,x,x2);z3=subs(F,x,x3); z4=subs(F,x,x4); b=z1;z2;z3;z4; g=ab; c=g(4,1);g(3,1);g(1,1
20、); y1=poly2sym(c)y1 =(4989443987306157*x2)/9007199254740992+(2546480093808581*x)/2251799813685248 + 4454695378303483/4503599627370496 %第一次迭代得到的多项式 X1 =-1 -1 -1 -1 ; % x1随迭代的变化X2 =-0.500000000000000 -0.438621271668792 -0.436958142528573 -0.436958064363172 ; % x2随迭代的变化X3 =0.500000000000000 0.560938664
21、934263 0.560058700967099 0.560057761761697 ; % x3随迭代的变化X4 =1 1 1 1 ; % x4随迭代的变化P1 = 0.367879441171442 0.412216302056878 0.412896544738883 0.412896829574159 ; % y1随迭代的变化P2 =0.606530659712633 0.562193798827198 0.599907881178971 0.600981082304536 ; %y2随迭代的变化P3 =1.648721270700128 1.693058131585564 1.797
22、333670243964 1.795792657883955 ; %y3随迭代的变化P4 =2.718281828459046 2.673944967573610 2.673264724891605 2.673264440056329 ; % y4随迭代的变化 x=-1:0.1:1;y1=(4989443987306157*x.2)/9007199254740992+(2546480093808581*x)/2251799813685248+4454695378303483/4503599627370496 ;y=1247589209708019/2251799813685248*x.2+50
23、89895364143919/4503599627370496*x+8908477905081045/9007199254740992;plot(x,exp(x),x,y1,+,x,y,-.)grid onlegend(原函数,初始误差曲线,最终误差曲线)hold onplot(X1,P1,v,X2,P2,v,X3,P3,v,X4,P4,v)图3第三题 最小二乘法为了制定生产计划，某羊毛衫厂记录了一部一年来羊毛衫的销售量，按月份得到表（1）的数据表，其销售量单位为箱，试建立月份（x）和销售（y）之间的关系。月份（x）1 2 3 4 5 6 7 8 9 10 11 12销量（y）256 201
24、159 61 77 40 17 25 103 156 222 345表1a=1 1 1;1 2 4;1 3 9; 1 4 16;1 5 25;1 6 36;1 7 49;1 8 64;1 9 81;1 10 100;1 11 121;1 12 144;b=256;201;159;61;77;40;17;25;103;156;222;345;x=inv(a*a)*a*bx = 386.0000 -113.42669.0420x1=1:0.01:12;y=9.0420.*x1.2-113.4266.*x1+386;plot(x1,y)hold onx2=1,2,3,4,5,6,7,8,9,10,1
25、1,12;y2=256,201,159,61,77,40,17,25,103,156,222,345;plot(x2,y2,.)图4第四题 用复化梯形公式计算积分。问积分区间要多少等分才能保证计算结果有5位有效数字？建立Trapezoid.m文件function I,step= Trapezoid (f,a,b,eps) %f为函数，a为积分上限，b为积分下限，eps为积分精度，step为划分区间个数if(nargin=3) eps=1.0e-4;endn=1;h=(b-a)/2;I1=2;I2=(subs(sym(f),findsym(sym(f),a)+subs(sym(f),findsym(sym(f),b)/h;while abs(I2-I1)eps n=n+1; h=(b-a)/n; I1=I2; I2=0; for i=0:n-1 %第n次的复化梯形公式积分 x=a+h*i; %i=0和n-1时，分别代表积分区间的左右端点 x1=x+h;I2=I2+(h/2)*(subs(sym(f),findsym(sym(f),x)+subs(sym(f),findsym(sym(f),x1); endendI=I2;step=n;输入:q,s=Trapezoid (exp(x),0,1,1e-6)q = 1.7183s = 67
