第六章函數(shù)的插值方法

1、第三篇 第六章 函數(shù)的插值方法第六章 函數(shù)的插值方法6.1 插值問(wèn)題及其誤差 6.1.2 與插值有關(guān)的MATLAB 函數(shù)(一) POLY2SYM 函數(shù)調(diào)用格式一：poly2sym (C)調(diào)用格式二：f1=poly2sym(C,'V') 或 f2=poly2sym(C, sym ('V') ),(二) POLYVAL 函數(shù)調(diào)用格式：Y = polyval(P,X)(三) POLY 函數(shù)調(diào)用格式：Y = poly (V)(四) CONV 函數(shù)調(diào)用格式：C =conv (A, B)例 6.1.2 求三個(gè)一次多項(xiàng)式、和的積.它們的零點(diǎn)分別依次為0.4,0.8,1.2.解

2、 我們可以用兩種MATLAB程序求之.方法1 如輸入MATLAB程序>> X1=0.4,0.8,1.2; l1=poly(X1), L1=poly2sym (l1)運(yùn)行后輸出結(jié)果為l1 = 1.0000 -2.4000 1.7600 -0.3840L1 = x3-12/5*x2+44/25*x-48/125方法2 如輸入MATLAB程序>> P1=poly(0.4);P2=poly(0.8);P3=poly(1.2); C =conv (conv (P1, P2), P3) , L1=poly2sym (C)運(yùn)行后輸出的結(jié)果與方法1相同.(五) DECONV 函數(shù)調(diào)用格

3、式：Q,R =deconv (B,A)(六) roots(poly(1:n)命令調(diào)用格式：roots(poly(1:n) (七) det(a*eye(size (A) - A)命令調(diào)用格式：b=det(a*eye(size (A) - A)6.2 拉格朗日（Lagrange）插值及其MATLAB程序6.2.1 線性插值及其MATLAB程序例6.2.1 已知函數(shù)在上具有二階連續(xù)導(dǎo)數(shù)，且滿足條件.求線性插值多項(xiàng)式和函數(shù)值，并估計(jì)其誤差.解 輸入程序>> X=1,3;Y=1,2; l01= poly(X(2)/( X(1)- X(2), l11= poly(X(1)/( X(2)- X(

4、1), l0=poly2sym (l01),l1=poly2sym (l11), P = l01* Y(1)+ l11* Y(2), L=poly2sym (P),x=1.5; Y = polyval(P,x)運(yùn)行后輸出基函數(shù)l0和l1及其插值多項(xiàng)式的系數(shù)向量P（略）、插值多項(xiàng)式L和插值Y為l0 = l1 = L = Y =-1/2*x+3/2 1/2*x-1/2 1/2*x+1/2 1.2500輸入程序>> M=5;R1=M*abs(x-X(1)* (x-X(2)/2運(yùn)行后輸出誤差限為 R1 = 1.8750例6.2.2 求函數(shù)e在上線性插值多項(xiàng)式，并估計(jì)其誤差.解 輸入程序&g

5、t;> X=0,1; Y =exp(-X) , 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),運(yùn)行后輸出基函數(shù)l0和l1及其插值多項(xiàng)式的系數(shù)向量P和插值多項(xiàng)式L為l0 = l1 = P =-x+1 x -0.6321 1.0000L =-1423408956596761/2251799813685248*x+1 輸入程序>> M=1;x=0:0.001:

6、1; R1=M*max(abs(x-X(1).*(x-X(2)./2運(yùn)行后輸出誤差限為 R1 = 0.1250.6.2.2 拋物線插值及其MATLAB程序例6.2.3 求將區(qū)間 0, /2 分成等份，用產(chǎn)生個(gè)節(jié)點(diǎn)，然后根據(jù)（6.9）和（6.13）式分別作線性插值函數(shù)和拋物線插值函數(shù).用它們分別計(jì)算cos (/6) (取四位有效數(shù)字)，并估計(jì)其誤差.解 輸入程序>> X=0,pi/2; Y =cos(X) ,l01= poly(X(2)/( X(1)- X(2), l11= poly(X(1)/( X(2)- X(1), l0=poly2sym (l01),l1=poly2sym (

7、l11), P = l01* Y(1)+ l11* Y(2), L=poly2sym (P),x=pi/6; Y = polyval(P,x)運(yùn)行后輸出基函數(shù)l0和l1及其插值多項(xiàng)式的系數(shù)向量P、插值多項(xiàng)式和插值為l0 =-5734161139222659/9007199254740992*x+1l1 =5734161139222659/9007199254740992*xP = -0.6366 1.0000L =-5734161139222659/9007199254740992*x+1Y =0.6667輸入程序>> M=1;x=pi/6; R1=M*abs(x-X(1)*(x-

8、X(2)/2運(yùn)行后輸出誤差限為R1 = 0.2742.（2） 輸入程序>> X=0:pi/4:pi/2; Y =cos(X) ,l01= conv (poly(X(2),poly(X(3)/( X(1)- X(2)* ( X(1)- X(3), l11= conv (poly(X(1), poly(X(3)/( X(2)- X(1)* ( X(2)- X(3),l21= conv (poly(X(1), poly(X(2)/( X(3)- X(1)* ( X(3)- X(2),l0=poly2sym (l01),l1=poly2sym (l11),l2=poly2sym (l21)

9、,P = l01* Y(1)+ l11* Y(2) + l21* Y(3), L=poly2sym (P),x=pi/6; Y = polyval(P,x)運(yùn)行后輸出基函數(shù)l01、l11和l21及其插值多項(xiàng)式的系數(shù)向量P、插值多項(xiàng)式L和插值Y為l0 =228155022448185/281474976710656*x2-2150310427208497/1125899906842624*x+1l1 =-228155022448185/140737488355328*x2+5734161139222659/2251799813685248*xl2 =228155022448185/2814749

10、76710656*x2-5734161139222659/9007199254740992*xP = -0.3357 -0.1092 1.0000L=-6048313895780875/18014398509481984*x2-7870612110600739/72057594037927936*x+1Y = 0.8508輸入程序>> M=1;x=pi/6; R2=M*abs(x-X(1)*(x-X(2) *(x-X(3)/6運(yùn)行后輸出誤差限為R2 =0.0239.6.2.3 次拉格朗日（Lagrange）插值及其MATLAB程序例6.2.4 給出節(jié)點(diǎn)數(shù)據(jù)，作三次拉格朗日插值多項(xiàng)式

11、計(jì)算，并估計(jì)其誤差.解 輸入程序>> X=-2,0,1,2; Y =17,1,2,17;p1=poly(X(1); p2=poly(X(2);p3=poly(X(3); p4=poly(X(4); l01= conv ( conv (p2, p3), p4)/( X(1)- X(2)* ( X(1)- X(3) * ( X(1)- X(4), l11= conv ( conv (p1, p3), p4)/( X(2)- X(1)* ( X(2)- X(3) * ( X(2)- X(4),l21= conv ( conv (p1, p2), p4)/( X(3)- X(1)* ( X

12、(3)- X(2) * ( X(3)- X(4),l31= conv ( conv (p1, p2), p3)/( X(4)- X(1)* ( X(4)- X(2) * ( X(4)- X(3),l0=poly2sym (l01),l1=poly2sym (l11),l2=poly2sym (l21), l3=poly2sym (l31),P = l01* Y(1)+ l11* Y(2) + l21* Y(3) + l31* Y(4),運(yùn)行后輸出基函數(shù)l0，l1，l2和l3及其插值多項(xiàng)式的系數(shù)向量P（略）為l0 =-1/24*x3+1/8*x2-1/12*x，l1 =1/4*x3-1/4*x2

13、-x+1l2 =-1/3*x3+4/3*x，l3 =1/8*x3+1/8*x2-1/4*x輸入程序>> L=poly2sym (P),x=0.6; Y = polyval(P,x)運(yùn)行后輸出插值多項(xiàng)式和插值為L(zhǎng) = Y =x3+4*x2-4*x+1 0.2560.輸入程序>> syms M; x=0.6; R3=M*abs(x-X(1)*(x-X(2) *(x-X(3) *(x-X(4)/24運(yùn)行后輸出誤差限為R3 =91/2500*M即 R3 ， .6.2.5 拉格朗日多項(xiàng)式和基函數(shù)的MATLAB程序求拉格朗日插值多項(xiàng)式和基函數(shù)的MATLAB主程序function C

14、, L,L1,l=lagran1(X,Y)m=length(X); L=ones(m,m);for k=1: m V=1; for i=1:m if k=i V=conv(V,poly(X(i)/(X(k)-X(i); endendL1(k,:)=V; l(k,:)=poly2sym (V)endC=Y*L1;L=Y*l例6.2.5 給出節(jié)點(diǎn)數(shù)據(jù)， ，作五次拉格朗日插值多項(xiàng)式和基函數(shù)，并寫(xiě)出估計(jì)其誤差的公式.解 在MATLAB工作窗口輸入程序>> X=-2.15 -1.00 0.01 1.02 2.03 3.25;Y=17.03 7.24 1.05 2.03 17.06 23.05

15、;C, L ,L1,l= lagran1(X,Y)運(yùn)行后輸出五次拉格朗日插值多項(xiàng)式L及其系數(shù)向量C，基函數(shù)l及其系數(shù)矩陣L1如下C = -0.2169 0.0648 2.1076 3.3960 -4.5745 1.0954L =1.0954-4.5745*x+3.3960*x2+2.1076*x3+0.0648*x4-0.2169*x5L1 = -0.0056 0.0299 -0.0323 -0.0292 0.0382 -0.0004 0.0331 -0.1377 -0.0503 0.6305 -0.4852 0.0048 -0.0693 0.2184 0.3961 -1.2116 -0.31

16、66 1.0033 0.0687 -0.1469 -0.5398 0.6528 0.9673 -0.0097 -0.0317 0.0358 0.2530 -0.0426 -0.2257 0.0023 0.0049 0.0004 -0.0266 0.0001 0.0220 -0.0002l = -0.0056*x5+0.0299*x4-0.0323*x3-0.0292*x2+0.0382*x-0.0004 0.0331*x5-0.1377*x4-0.0503*x3+0.6305*x2-0.4852*x+0.0048 -0.0693*x5+0.2184*x4+0.3961*x3-1.2116*x2

17、-0.3166*x+1.0033 0.0687*x5-0.1469*x4-0.5398*x3+0.6528*x2+0.9673*x-0.0097 -0.0317*x5+0.0358*x4+0.2530*x3-0.0426*x2-0.2257*x+0.0023 0.0049*x5+0.0004 *x4-0.0266*x3+0.0001*x2+0.0220*x-0.0002估計(jì)其誤差的公式為，.6.2.6 拉格朗日插值及其誤差估計(jì)的MATLAB程序拉格朗日插值及其誤差估計(jì)的MATLAB主程序function y,R=lagranzi(X,Y,x,M)n=length(X); m=length(x)

18、;for i=1:m z=x(i);s=0.0; for k=1:n p=1.0; q1=1.0; c1=1.0;for j=1:n if j=kp=p*(z-X(j)/(X(k)-X(j); end q1=abs(q1*(z-X(j);c1=c1*j; end s=p*Y(k)+s; end y(i)=s;endR=M*q1/c1;例 6.2.6 已知，用拉格朗日插值及其誤差估計(jì)的MATLAB主程序求的近似值，并估計(jì)其誤差.解 在MATLAB工作窗口輸入程序>> x=2*pi/9; M=1; X=pi/6 ,pi/4, pi/3;Y=0.5,0.7071,0.8660; y,R=

19、lagranzi(X,Y,x,M)運(yùn)行后輸出插值y及其誤差限R為y = R =0.6434 8.8610e-004.6.3 牛頓（Newton）插值及其MATLAB程序6.3.3 牛頓插值多項(xiàng)式、差商和誤差公式的MATLAB程序求牛頓插值多項(xiàng)式和差商的MATLAB主程序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-

20、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;例6.3.3 給出節(jié)點(diǎn)數(shù)據(jù)，,作五階牛頓插值多項(xiàng)式和差商，并寫(xiě)出其估計(jì)誤差的公式.解 （1）保存名為newpoly.m的M文件.（2）輸入MATL

21、AB程序>> X=-2.15 -1.00 0.01 1.02 2.03 3.25;Y=17.03 7.24 1.05 2.03 17.06 23.05; A,C,L,wcgs,Cw= newdcw (X,Y)運(yùn)行后輸出差商矩陣A,五階牛頓插值多項(xiàng)式L及其系數(shù)向量C, 插值余項(xiàng)公式L及其向量Cw如下A = 17.0300 0 0 0 0 0 7.2400 -8.5130 0 0 0 0 1.0500 -6.1287 1.1039 0 0 0 2.0300 0.9703 3.5144 0.7604 0 0 17.0600 14.8812 6.8866 1.1129 0.0843 0 2

22、3.0500 4.9098 -4.4715 -3.5056 -1.0867 -0.2169C = -0.2169 0.0648 2.1076 3.3960 -4.5745 1.0954L =-7813219284746629/36028797018963968*x5+583849564517807/9007199254740992*x4+593245028711263/281474976710656*x3+3823593773002357/1125899906842624*x2-321902673270315/70368744177664*x+308328649211299/281474976

23、710656wcgs =1/720*M*(x6-79/25*x5-14201/2500*x4+4934097026900981/281474976710656*x3+154500712237335/35184372088832*x2-8170642380559269/562949953421312*x+5212760744134241/36028797018963968)Cw = 0.0014 -0.0044 -0.0079 0.0243 0.0061 -0.0202 0.0002即L =1.0954-4.5745*x+3.3960*x2+2.1076*x3+0.0648*x4-0.2169*

24、x5.估計(jì)其誤差的公式為，.例6.3.4 求函數(shù)e在上五階牛頓插值多項(xiàng)式，估計(jì)其誤差的公式和誤差限公式.用它們計(jì)算，并估計(jì)其誤差.解 （1）輸入MATLAB程序>> X=2:4/5:6; Y=-7*exp(-X/5); A,C,L,wcgs,Cw= newploy(X,Y), x1=2:0.001:6; M=max(-7*exp(-x1/5)/(56),運(yùn)行后輸出差商矩陣A, 五階牛頓插值多項(xiàng)式L及其系數(shù)向量C, 插值余項(xiàng)公式L及其向量Cw如下A = -4.6922 0 0 0 0 0 -3.9985 0.8672 0 0 0 0 -3.4073 0.7390 -0.0801 0

25、0 0 -2.9035 0.6297 -0.0683 0.0049 0 0 -2.4742 0.5366 -0.0582 0.0042 -0.0002 0 -2.1084 0.4573 -0.0496 0.0036 -0.0002 0.0000C = 0.0000 -0.0004 0.0089 -0.1389 1.3985 -6.9991L =9721799720875/1152921504606846976*x5-3503994098647815/9223372036854775808*x4+160742008798419/18014398509481984*x3-1251152213853

26、501/9007199254740992*x2+6298131904328647/4503599627370496*x-3940156929554013/562949953421312wcgs =1/720*M*(x6-24*x5+1172/5*x4-5952/5*x3+7276634802928539/2199023255552*x2-5237461186650519/1099511627776*x+6085939356447121/2199023255552)Cw =0.0014 -0.0333 0.3256 -1.6533 4.5959 -6.6159 3.8438M = -1.3494

27、e-004（2）輸入MATLAB程序>> syms xwcgs1=1/720*M*1/720*M*(x6-24*x5+1172/5*x4-5952/5*x3+7276634802928539/2199023255552*x2-5237461186650519/1099511627776*x+6085939356447121/2199023255552),運(yùn)行后輸出誤差限公式wcgs1如下wcgs1 =5565367633581443/158456325028528675187087900672*x6-16696102900744329/19807040628566084398385

28、987584*x5+1630652716639362799/198070406285660843983859875840*x4-517579189923074199/12379400392853802748991242240*x3+40497147813610772932297365501777/348449143727040986586495598010130648530944*x2-29148396970323855270001164718917/174224571863520493293247799005065324265472*x+338704899143102839266652763

29、75603/348449143727040986586495598010130648530944（3）輸入MATLAB程序>> x=3.156;y=9721799720875/1152921504606846976*x5-3503994098647815/9223372036854775808*x4+160742008798419/18014398509481984*x3-1251152213853501/9007199254740992*x2+6298131904328647/4503599627370496*x-3940156929554013/562949953421312w

30、cgs2=1/720*M*(x6-24*x5+1172/5*x4-5952/5*x3+7276634802928539/2199023255552*x2-5237461186650519/1099511627776*x+6085939356447121/2199023255552)運(yùn)行后輸出的近似值y，及其誤差限wcgs2如下y = -3.7237wcgs2 = -2.4764e-0076.3.4 牛頓插值及其誤差估計(jì)的MATLAB程序牛頓插值及其誤差估計(jì)的MATLAB主程序function y,R= newcz(X,Y,x,M)n=length(X); m=length(x);for t=1

31、:m z=x(t); A=zeros(n,n);A(:,1)=Y' s=0.0; p=1.0; q1=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 q1=abs(q1*(z-X(j-1);c1=c1*j; end C=A(n,n);q1=abs(q1*(z-X(n);for k=(n-1):-1:1C=conv(C,poly(X(k);d=length(C); C(d)=C(d)+A(k,k);end y(k)= polyval(C, z);endR=M*q1/c1;例

32、6.3.5 已知，用牛頓插值法求的近似值，估計(jì)其誤差，并與例 6.2.6的計(jì)算結(jié)果比較.解 方法1（牛頓插值及其誤差估計(jì)的MATLAB主程序）輸入MATLAB程序>> x=2*pi/9;M=1; X=pi/6 ,pi/4, pi/3; Y=0.5,0.7071,0.8660; y,R= newcz(X,Y,x,M)運(yùn)行后輸出y = R =0.6434 8.8610e-004方法2 （求牛頓插值多項(xiàng)式和差商的MATLAB主程序）輸入MATLAB程序>> x=2*pi/9; X=pi/6 ,pi/4, pi/3; Y=0.5,0.7071,0.8660; M=1; A,C,

33、L,wcgs,Cw= newploy(X,Y), y=polyval(C,x)運(yùn)行后輸出結(jié)果A = 0.5000 0 0 0.7071 0.7911 0 0.8660 0.6070 -0.3516C = -0.3516 1.2513 -0.0588L =-1583578379808357/4503599627370496*x2+1408883391907005/1125899906842624*x-132405829044691/2251799813685248wcgs =1/6*M*(x3-3/4*x2*pi+4012734077357799/2251799813685248*x-77577

34、69783530263/18014398509481984)Cw = 0.1667 -0.3927 0.2970 -0.0718y = 0.6434上述兩種方法計(jì)算y的結(jié)果相同.6.3.5 牛頓插值法的MATLAB綜合程序求牛頓插值多項(xiàng)式、差商、插值及其誤差估計(jì)的MATLAB主程序function y,R,A,C,L=newdscg(X,Y,x,M)n=length(X); m=length(x);for t=1:m z=x(t); A=zeros(n,n);A(:,1)=Y's=0.0; p=1.0; q1=1.0; c1=1.0; for j=2:n for i=j:n A(i,j

35、)=(A(i,j-1)- A(i-1,j-1)/(X(i)-X(i-j+1); end q1=abs(q1*(z-X(j-1);c1=c1*j; end C=A(n,n);q1=abs(q1*(z-X(n);for k=(n-1):-1:1C=conv(C,poly(X(k);d=length(C);C(d)=C(d)+A(k,k);end y(k)= polyval(C, z);endR=M*q1/c1;L(k,:)=poly2sym(C);例6.3.6 給出節(jié)點(diǎn)數(shù)據(jù)，作三階牛頓插值多項(xiàng)式，計(jì)算，并估計(jì)其誤差.解 首先將名為newdscg.m的程序保存為M文件，然后在MATLAB工作窗口輸入

36、程序>> syms M,X=-4,0,1,2; Y =27,1,2,17; x=-2.345; y,R,A,C,P=newdscg(X,Y,x,M)運(yùn)行后輸出插值y及其誤差限公式R，三階牛頓插值多項(xiàng)式P及其系數(shù)向量C，差商的矩陣A如下y = 22.3211R =1323077530165133/562949953421312*M（即R =2.3503*M）A= 27.0000 0 0 0 1.0000 -6.5000 0 0 2.0000 1.0000 1.5000 0 17.0000 15.0000 7.0000 0.9167C = 0.9167 4.2500 -4.1667 1

37、.0000P =11/12*x3+17/4*x2-25/6*x+1例6.3.7 將區(qū)間 0,/2 分成等份，用產(chǎn)生個(gè)節(jié)點(diǎn)，求二階和三階牛頓插值多項(xiàng)式和.用它們分別計(jì)算sin (/7) (取四位有效數(shù)字)，并估計(jì)其誤差.解 首先將名為newdscg.m的程序保存為M文件，然后在MATLAB工作窗口輸入程序>> X1=0:pi/4:pi/2; Y1 =sin(X1); M=1; x=pi/7; X2=0:pi/6:pi/2; Y2 =sin(X2); y1,R1,A1,C1,P2=newdscg(X1,Y1,x,M), y2,R2,A2,C2,P3=newdscg(X2,Y2,x,M)

38、運(yùn)行后輸出插值y1=和y2=及其誤差限R1和R2，二階和三階牛頓插值多項(xiàng)式P2和P3及其系數(shù)向量C1和C2，差商的矩陣A1和A2如下y1 = 0.4548R1 = 0.0282A1 = 0 0 0 0.7071 0.9003 0 1.0000 0.3729 -0.3357C1 = -0.3357 1.1640 0P2 =-3024156947890437/9007199254740992*x2+163820246322191/140737488355328*xy2 = 0.4345R2 = 9.3913e-004A2 = 0 0 0 0 0.5000 0.9549 0 0 0.8660 0.6

39、991 -0.2443 0 1.0000 0.2559 -0.4232 -0.1139C2 = -0.1139 -0.0655 1.0204 0P3 =-1025666884451963/9007199254740992*x3-4717668559161127/72057594037927936*x2+4595602396951547/4503599627370496*x6.4 埃爾米特（Hermite）插值及其MATLAB程序6.4.3 埃爾米特插值多項(xiàng)式和誤差公式的MATLAB程序求埃爾米特插值多項(xiàng)式和誤差公式的MATLAB主程序function Hc, Hk,wcgs,Cw= hermi

40、te (X,Y,Y1)m=length(X); n=M1;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+1if k=is=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);endH=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=

41、poly2sym(q);endfor i=1:t c1=c1*i; endsyms M,wcgs=M*Q/c1; Cw=q/c1;例6.4.3 給定函數(shù)在點(diǎn)處的函數(shù)值, ，和導(dǎo)數(shù)值，且，求函數(shù)在點(diǎn)處的五階埃爾米特插值多項(xiàng)式和誤差公式，計(jì)算并估計(jì)其誤差.解 （1）保存名為hermite.m的M文件.（2）在MATLAB工作窗口輸入程序>>X=pi/6,pi/4,pi/2; Y=0.5,0.7071,1;Y1=0.8660,0.7071,0; Hc, Hk,wcgs,Cw= hermite (X,Y,Y1)運(yùn)行后輸出五階埃爾米特插值多項(xiàng)式Hk及其系數(shù)向量Hc，誤差公式wcgs及其系數(shù)向

42、量Cw如下Hc = 1.0e+003 * 0.1912 -0.9273 1.6903 -1.4380 0.5751 -0.0866Hk =6725828781679091/35184372088832*x5-4078286086775209/4398046511104*x4+7434035571017927/4398046511104*x3-3162205449085973/2199023255552*x2+5058863928652835/8796093022208*x-6094057839958843/70368744177664wcgs =1/720*M*(x6-11/6*x5*pi+7

43、446708432019761/562949953421312*x4-4363745503235773/281474976710656*x3+21569239021155/2199023255552*x2-7178073637328281/2251799813685248*x+3758430567659515/9007199254740992)Cw =0.0014 -0.0080 0.0184 -0.0215 0.0136 -0.0044 0.0006當(dāng)時(shí)的誤差公式為R=0.001 4* x6-0.008 0*x5+0.018 4*x4-0.021 5*x3+0.013 6*x2-0.004

44、4*x+0.000 6（3）在MATLAB工作窗口輸入程序>> x=1.567;M=1; Hk=6725828781679091/35184372088832*x5-4078286086775209/4398046511104*x4+7434035571017927/4398046511104*x3-3162205449085973/2199023255552*x2+5058863928652835/8796093022208*x-6094057839958843/70368744177664,wcgs=1/720*M*(x6-11/6*x5*pi+7446708432019761

45、/562949953421312*x4-4363745503235773/281474976710656*x3+21569239021155/2199023255552*x2-7178073637328281/2251799813685248*x+3758430567659515/9007199254740992)運(yùn)行后輸出的近似值Hk及其誤差wcgs如下Hk = 2.5265wcgs = 1.3313e-0086.5 分段插值及其MATLAB程序6.5.1 高次插值的振蕩和MATLAB程序例6.5.1 作下列函數(shù)在指定區(qū)間上的次拉格朗日插值多項(xiàng)式 的圖形，并討論插值多項(xiàng)式的次數(shù)與誤差的關(guān)系.

46、（1），；（2），.解 將計(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（1）現(xiàn)提供作次拉格朗日插值多項(xiàng)式的圖形的MATLAB程序，保存名為 Li1.m 的M文件m=150; x=-pi:2*pi/(m-1): pi; y=tan(

47、cos(3(1/2)+sin(2*x)./(3+4*x.2);plot(x,y,'k-'),gtext('y=tan(cos(sqrt(3)+sin(2x)/(3+4x2)'),pausen=3; x0=-pi:2*pi/(3-1):pi; y0=tan(cos(3(1/2)+sin(2*x0)./(3+4*x0.2);y1=lagr1(x0,y0,x);hold on,plot(x,y1,'g<'),gtext('n=2'),pause,hold offn=5; x0=-pi:2*pi/(n-1):pi; y0=tan(c

48、os(3(1/2)+sin(2*x0)./(3+4*x0.2);y2=lagr1(x0,y0,x);hold on,plot(x,y2,'b:'),gtext('n=4'),pause,hold offn=7; x0=-pi:2*pi/(n-1):pi; y0=tan(cos(3(1/2)+sin(2*x0)./(3+4*x0.2);y3=lagr1(x0,y0,x);hold on,plot(x,y3,'rp'),gtext('n=6'),pause,hold offn=9; x0=-pi:2*pi/(n-1):pi; y0=t

49、an(cos(3(1/2)+sin(2*x0)./(3+4*x0.2);y4=lagr1(x0,y0,x);hold on,plot(x,y4,'m*'),gtext('n=8'),pause,hold offn=11; x0=-pi:2*pi/(n-1):pi; y0=tan(cos(3(1/2)+sin(2*x0)./(3+4*x0.2);y5=lagr1(x0,y0,x);hold on,plot(x,y5,'g:'),gtext('n=10')title('高次拉格朗日插值的振蕩')在MATLAB工作窗口輸

50、入名為 Li1.m的M文件的文件名 >> Li1.m回車(chē)運(yùn)行后，便會(huì)逐次畫(huà)出在區(qū)間上的次拉格朗日插值多項(xiàng)式 的圖形(略).（2）在MATLAB工作窗口輸入程序m=101; x=-5:10/(m-1):5; y=1./(1+x.2); z=0*x;plot(x,z,'r',x,y,'k-'),gtext('y=1/(1+x2)'),pausen=3; x0=-5:10/(n-1):5; y0=1./(1+x0.2); y1=lagr1(x0,y0,x);hold on,plot(x,y1,'g'),gtext('

51、n=2'),pause,hold offn=5; x0=-5:10/(n-1):5; y0=1./(1+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+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./

52、(1+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+x0.2); y5=lagr1(x0,y0,x);hold on,plot(x,y5,'m'),gtext('n=10')title('高次拉格朗日插值的振蕩')回車(chē)運(yùn)行后，便會(huì)逐次畫(huà)出在區(qū)間上的次拉格朗日插值多項(xiàng)式 的圖形(略). 6.5.4 作有關(guān)分段線性插值圖形的MATLAB程序作有關(guān)分段線性插值圖形的MATLAB主程序function s= xxczhjt1(x0,y0,xi,x,y)s= interp1(x0,