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

楊穎 應(yīng)用數(shù)學131802班 2013180502286.2 1.已知函數(shù)在1,7上具有二階連續(xù)導數(shù)，且滿足條件f(x)=1,f(7)=12,.求線性插值多項式和函數(shù)值f(3.5)，并估計其誤差.解： 輸入程序 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)運行后輸出基函數(shù)l0和l1及其插值多項式的系數(shù)向量P（略）、插值多項式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é)點數(shù)據(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)運行后輸出差商矩陣A,五階牛頓插值多項式L及其系數(shù)向量C, 插值余項公式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估計其誤差的公式為=（x+3.15)(x+1.00)(x-0.01)(x-1.02)(x-2.03)(x-3.25)，6.41.給定函數(shù)在點x=,處的函數(shù)值, ，和導數(shù)值，且，求函數(shù)在點處的三階埃爾米特插值多項式和誤差公式。解 （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)運行后輸出三階埃爾米特插值多項式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當=M時的誤差公式為R=0.0417-0.1091x+ 0.1057x-0.0449x+0.00706.51. 作函數(shù)在區(qū)間-5,5上的次拉格朗日插值多項式 的圖形，并討論插值多項式的次數(shù)與誤差的關(guān)系.解 將計算次拉格朗日插值多項式的值的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(高次拉格朗日插值的振蕩)回車運行后，便會逐次畫出在區(qū)間上的次拉格朗日插值多項式 的圖形從圖中可以看出，對于較大的|x|,隨著n的增大，L(x)振蕩越來越大，事實上可以證明，僅當是某一個正數(shù)）時，才有=g(x).而在此區(qū)間外，是發(fā)散的。由于在大范圍內(nèi)使用高次插值多項式逼近的效果往往不理想，從而促使人們轉(zhuǎn)而尋求簡單的低次插值方法。從圖中可以看出，對于較大的|x|,隨著n的增大，L(x)振蕩越來越大，事實上可以證明，僅當是某一個正數(shù)）時，才有=g(x).而在此區(qū)間外，是發(fā)散的。由于在大范圍內(nèi)使用高次插值多項式逼近的效果往往不理想，從而促使人們轉(zhuǎn)而尋求簡單的低次插值方法。6.7設(shè)函數(shù)定義在區(qū)間上，取n=10，按等距節(jié)點構(gòu)造四種插值函數(shù)S(x).(1) 用MATLAB程序計算各小區(qū)間中點處，作出及其后三種插值函數(shù)，作出節(jié)點、插值點、和四種插值函數(shù)的圖形；(2) 并用MATLAB程序計算各小區(qū)間中點處四種插值及其相對誤差；(3) 用MATLAB程序估計和S(x)在區(qū)間上的的誤差限。解 （1）記節(jié)點的橫坐標,插值點,i= 編寫并保存名為sanli679.m的M文件如下h記節(jié)點的橫坐標插值點，.= 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é)點,三次樣條插值點,被插值函數(shù),拉格朗日插值函數(shù),三次樣條函數(shù),分段三次埃爾米特插值函數(shù))title(y=tan(cos(sqrt(3)+sin(2x)/(3+4x2)及其三種插值函數(shù),節(jié)點和插值點的圖形)運行程序 n=10;sanli679得到各小區(qū)間中點處的函數(shù)值，分段線性插值、拉格朗日插值、三次樣條插值和分段埃爾米特插值及其絕對誤差、相對誤差和平均值的結(jié)果,作出及其后三種插值函數(shù)，插值點和節(jié)點的圖形(略). 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