數(shù)值分析,任玉杰教程課后題.doc

楊穎 應(yīng)用數(shù)學(xué)131802班 2013180502286.2 1.已知函數(shù)在1,7上具有二階連續(xù)導(dǎo)數(shù)，且滿足條件f(x)=1,f(7)=12,.求線性插值多項(xiàng)式和函數(shù)值f(3.5)，并估計(jì)其誤差.解： 輸入程序 X=1,7;Y=1,12; l01= poly(X(2)/( X(1)- X(2), l11= poly(X(1)/( X(2)- X(1), l0=poly2sym (l01),l1=poly2sym (l11), P = l01* Y(1)+ l11* Y(2), L=poly2sym (P),x=3.5; Y = polyval(P,x)運(yùn)行后輸出基函數(shù)l0和l1及其插值多項(xiàng)式的系數(shù)向量P（略）、插值多項(xiàng)式L和插值Y為l01 = -0.1667 1.1667l11 = 0.1667 -0.1667l0 =7/6 - x/6l1 =x/6 - 1/6P =1.8333 -0.8333L =(11*x)/6 - 5/6Y =5.58336.35.給出節(jié)點(diǎn)數(shù)據(jù) ，,作五階牛頓插值多項(xiàng)式和差商，并寫出其估計(jì)誤差的公式.解 （1）保存名為newpoly.m的M文件.function A,C,L,wcgs,Cw= newploy(X,Y)n=length(X); A=zeros(n,n); A(:,1)=Y; s=0.0; p=1.0; q=1.0; c1=1.0; for j=2:n for i=j:n A(i,j)=(A(i,j-1)- A(i-1,j-1)/(X(i)-X(i-j+1); end b=poly(X(j-1);q1=conv(q,b); c1=c1*j; q=q1; end C=A(n,n); b=poly(X(n); q1=conv(q1,b); for k=(n-1):-1:1C=conv(C,poly(X(k); d=length(C); C(d)=C(d)+A(k,k);endL(k,:)=poly2sym(C); Q=poly2sym(q1);syms Mwcgs=M*Q/c1; Cw=q1/c1;（2）輸入MATLAB程序 X=-3.15 -1.00 0.01 1.02 2.03 3.25;Y=37.03 7.24 1.05 2.03 17.06 23.05; A,C,L,wcgs,Cw= newpoly (X,Y)運(yùn)行后輸出差商矩陣A,五階牛頓插值多項(xiàng)式L及其系數(shù)向量C, 插值余項(xiàng)公式L及其向量Cw如下A = 37.0300 0 0 0 0 0 7.2400 -13.8558 0 0 0 0 1.0500 -6.1287 2.4453 0 0 0 2.0300 0.9703 3.5144 0.2564 0 0 17.0600 14.8812 6.8866 1.1129 0.1654 0 23.0500 4.9098 -4.4715 -3.5056 -1.0867 -0.1956C = -0.1956 -0.0479 2.2294 3.5063 -4.7185 1.0968L=-(3524262145772209*x5)/18014398509481984 - (431283038958129*x4)/9007199254740992 + (5020108164096803*x3)/2251799813685248 + (246736293371003*x2)/70368744177664 - (5312509994841131*x)/1125899906842624 + 154365339321353/140737488355328wcgs =(M*(x6 - (54*x5)/25 - (6869*x4)/625 + (204645801275343*x3)/8796093022208 + (5397560409786347*x2)/562949953421312 - (5994280682419117*x)/281474976710656 + 7637300625126911/36028797018963968)/720Cw = 0.0014 -0.0030 -0.0153 0.0323 0.0133 -0.0296 0.0003即L=1.09683171-4.71845673x+3.50633362x+2.22937587x-0.04788204x-0.19563585x估計(jì)其誤差的公式為=（x+3.15)(x+1.00)(x-0.01)(x-1.02)(x-2.03)(x-3.25)，6.41.給定函數(shù)在點(diǎn)x=,處的函數(shù)值, ，和導(dǎo)數(shù)值，且，求函數(shù)在點(diǎn)處的三階埃爾米特插值多項(xiàng)式和誤差公式。解 （1）保存名為hermite.m的M文件.function Hc, Hk,wcgs,Cw= hermite (X,Y,Y1)m=length(X); n=1;s=0; H=0;q=1;c1=1; L=ones(m,m); G=ones(1,m);for k=1:n+1 V=1; for i=1:n+1 if k=i s=s+(1/(X(k)-X(i); V=conv(V,poly(X(i)/(X(k)-X(i); end h=poly(X(k); g=(1-2*h*s); G=g*Y(k)+h*Y1(k); end H=H+conv(G,conv(V,V); b=poly(X(k);b2=conv(b,b); q=conv(q,b2); t=2*n+2; Hc=H;Hk=poly2sym (H); Q=poly2sym(q);endfor i=1:t c1=c1*i; endsyms M,wcgs=M*Q/c1; Cw=q/c1;（2）在MATLAB工作窗口輸入程序X=pi/6,pi/4; Y=0.5,0.7071;Y1=0.8660,0.7071; Hc, Hk,wcgs,Cw= hermite (X,Y,Y1)運(yùn)行后輸出三階埃爾米特插值多項(xiàng)式Hk及其系數(shù)向量Hc，誤差公式wcgs及其系數(shù)向量Cw如下Hc =96.2945 -166.7424 94.7792 -16.9739Hk=(3388060314322531*x3)/35184372088832 - (5866727782923391*x2)/35184372088832 + (6669496320254561*x)/70368744177664 - 2388870952879233/140737488355328wcgs =(M*(x4 - (5*pi*x3)/6 + (1427607315983063*x2)/562949953421312 - (4848606114706415*x)/4503599627370496 + 6092938140044889/36028797018963968)/24Cw =0.0417 -0.1091 0.1057 -0.0449 0.0070當(dāng)=M時(shí)的誤差公式為R=0.0417-0.1091x+ 0.1057x-0.0449x+0.00706.51. 作函數(shù)在區(qū)間-5,5上的次拉格朗日插值多項(xiàng)式 的圖形，并討論插值多項(xiàng)式的次數(shù)與誤差的關(guān)系.解 將計(jì)算次拉格朗日插值多項(xiàng)式的值的MATLAB程序保存名為lagr1.m的M文件.function y=lagr1(x0,y0,x)n=length(x0); m=length(x);for i=1:m z=x(i);s=0.0;for k=1:n p=1.0;for j=1:nif j=k p=p*(z-x0(j)/(x0(k)-x0(j); end end s=p*y0(k)+s; end y(i)=s;End在MATLAB工作窗口輸入程序m=101; x=-5:10/(m-1):5; y=1./(1+25*x.2); z=0*x;plot(x,z,r,x,y,k-),gtext(y=1/(1+25*x2),pausen=3; x0=-5:10/(n-1):5; y0=1./(1+25*x0.2); y1=lagr1(x0,y0,x);hold on,plot(x,y1,g),gtext(n=2),pause,hold offn=5; x0=-5:10/(n-1):5; y0=1./(1+25*x0.2); y2=lagr1(x0,y0,x);hold on,plot(x,y2,b:),gtext(n=4),pause,hold offn=7; x0=-5:10/(n-1):5; y0=1./(1+25*x0.2); y3=lagr1(x0,y0,x);hold on,plot(x,y3,r),gtext(n=6),pause,hold offn=9; x0=-5:10/(n-1):5; y0=1./(1+25*x0.2); y4=lagr1(x0,y0,x);hold on,plot(x,y4,r:),gtext(n=8),pause,hold offn=11; x0=-5:10/(n-1):5; y0=1./(1+25*x0.2); y5=lagr1(x0,y0,x);hold on,plot(x,y5,m),gtext(n=10)title(高次拉格朗日插值的振蕩)回車運(yùn)行后，便會(huì)逐次畫出在區(qū)間上的次拉格朗日插值多項(xiàng)式 的圖形從圖中可以看出，對(duì)于較大的|x|,隨著n的增大，L(x)振蕩越來越大，事實(shí)上可以證明，僅當(dāng)是某一個(gè)正數(shù)）時(shí)，才有=g(x).而在此區(qū)間外，是發(fā)散的。由于在大范圍內(nèi)使用高次插值多項(xiàng)式逼近的效果往往不理想，從而促使人們轉(zhuǎn)而尋求簡(jiǎn)單的低次插值方法。從圖中可以看出，對(duì)于較大的|x|,隨著n的增大，L(x)振蕩越來越大，事實(shí)上可以證明，僅當(dāng)是某一個(gè)正數(shù)）時(shí)，才有=g(x).而在此區(qū)間外，是發(fā)散的。由于在大范圍內(nèi)使用高次插值多項(xiàng)式逼近的效果往往不理想，從而促使人們轉(zhuǎn)而尋求簡(jiǎn)單的低次插值方法。6.7設(shè)函數(shù)定義在區(qū)間上，取n=10，按等距節(jié)點(diǎn)構(gòu)造四種插值函數(shù)S(x).(1) 用MATLAB程序計(jì)算各小區(qū)間中點(diǎn)處，作出及其后三種插值函數(shù)，作出節(jié)點(diǎn)、插值點(diǎn)、和四種插值函數(shù)的圖形；(2) 并用MATLAB程序計(jì)算各小區(qū)間中點(diǎn)處四種插值及其相對(duì)誤差；(3) 用MATLAB程序估計(jì)和S(x)在區(qū)間上的的誤差限。解 （1）記節(jié)點(diǎn)的橫坐標(biāo),插值點(diǎn),i= 編寫并保存名為sanli679.m的M文件如下h記節(jié)點(diǎn)的橫坐標(biāo)插值點(diǎn)，.= tan(cos(3(1/2)+sin(2*X)./(3+4*X.2);YL=lagr1(x0,y0,X); YS=interp1(x0,y0,X,spline);YH=interp1(x0,y0,X, pchip); yL=lagr1(x0,y0,x); yX=interp1(x0,y0,x); yS=interp1(x0,y0,x,spline); yH=interp1(x0,y0,x, pchip); RL=abs(y-yL)./y);RS=abs(y-yS)./y); RH=abs(y-yH)./y); RX=abs(y-yX)./y);RLj=abs(y-yL);mRLj=mean(RLj);RSj=abs(y-yS);mRSj=mean(RSj);RHj=abs(y-yH);RXj=abs(y-yX);mRHj=mean(RHj);mRXj=mean(RXj);mRL=mean(RL);mRX=mean(RX);mRS=mean(RS);mRH=mean(RH);CZ=x y yL yX yS yH,R=x RL RX RS RH,mR=mRL mRX mRS mRHRj=x RLj RXj RSj RHj,mRj=mRLj mRXj mRSj mRHj,plot(x0,y0,bo,x,yS,r*,X,Y,k-,X,YL,g-.,X,YS,c:, X,YH,m-),legend(節(jié)點(diǎn),三次樣條插值點(diǎn),被插值函數(shù),拉格朗日插值函數(shù),三次樣條函數(shù),分段三次埃爾米特插值函數(shù))title(y=tan(cos(sqrt(3)+sin(2x)/(3+4x2)及其三種插值函數(shù),節(jié)點(diǎn)和插值點(diǎn)的圖形)運(yùn)行程序 n=10;sanli679得到各小區(qū)間中點(diǎn)處的函數(shù)值，分段線性插值、拉格朗日插值、三次樣條插值和分段埃爾米特插值及其絕對(duì)誤差、相對(duì)誤差和平均值的結(jié)果,作出及其后三種插值函數(shù)，插值點(diǎn)和節(jié)點(diǎn)的圖形(略). CZ = -2.8274 1.5499 0.9343 1.5483 1.5485 1.5487 -2.1991 1.5330 1.6668 1.5344 1.5350 1.5327 -1.5708 1.5269 1.4666 1.5293 1.5204 1.5293 -0.9425 1.5334 1.5848 1.5204 1.5585 1.5258 -0.3142 1.3800 1.3133 1.3099 1.3173 1.3070 0.3142 0.9803 1.0247 1.1057 1.0331 1.1032 0.9425 1.3060 1.2818 1.2785 1.2788 1.2600 1.5708 1.5269 1.5433 1.5032 1.5359 1.5202 2.1991 1.5553 1.5466 1.5530 1.5533 1.5544 2.8274 1.5556 1.4976 1.5553 1.5570 1.5559R = -2.8274 0.3972 0.0010 0.0009 0.0008 -2.1991 0.0873 0.0009 0.0013 0.0002 -1.5708 0.0394 0.0016 0.0042 0.0016 -0.9425 0.0335 0.0085 0.0164 0.0049 -0.3142 0.0483 0.0508 0.0454 0.0529 0.3142 0.0452 0.1279 0.0538 0.1253 0.9425 0.0186 0.0211 0.0209 0.0352 1.5708 0.0108 0.0155 0.0059 0.0043 2.1991 0.0056 0.0015 0.0013 0.0006 2.8274 0.0373 0.0002 0.0009 0.0002mR = 0.0723 0.0229 0.0151 0.0226Rj = -2.8274 0.6156 0.0016 0.0014 0.0012 -2.1991 0.1338 0.0014 0.0020 0.0003 -1.5708 0.0602 0.0025 0.0065 0.0025 -0.9425 0.0514 0.0130 0.0251 0.0076 -0.3142 0.0667 0.0701 0.0627 0.0730 0.3142 0.0443 0.1254 0.0528 0.1228 0.9425 0.0243 0.0275 0.0273 0.0460 1.5708 0.0165 0.0236 0.0090 0.0066 2.1991 0.0087 0.0023 0.0020 0.0009 2.8274 0.0580 0.0002 0.0014 0.0003mRj = 0.1080 0.0268 0.0190 0.0261(3).在MATLAB工作窗口輸入程序 syms x y=tan(cos(3(1/2)+sin(2*x)/(3+4*x2);yxx=dif