2022年第六章函数的插值方法

《2022年第六章函数的插值方法》由会员分享，可在线阅读，更多相关《2022年第六章函数的插值方法（29页珍藏版）》请在金锄头文库上搜索。

1、学习必备欢迎下载6.1 插值问题及其误差6.1.2 与插值有关的 MATLAB 函数( 一) POLY2SYM函数调用格式一 ：poly2sym (C) 调用格式二 ：f1=poly2sym(C,V) 或 f2=poly2sym(C, sym (V) ),( 二) POLYVAL函数调用格式 ：Y = polyval(P,X)( 三) POLY函数调用格式 ：Y = poly (V)( 四) CONV函数调用格式 ：C =conv (A, B)例 6.1.2求三个一次多项式)(xg、)(xh和)(xf的积)()()(xhxgxf.它们的零点分别依次为0.4,0.8,1.2. 解我们可以用两种M

2、ATLAB 程序求之 . 方法 1如输入 MATLAB 程序 X1=0.4,0.8,1.2; l1=poly(X1), L1=poly2sym (l1) 运行后输出结果为l1 = 1.0000 -2.4000 1.7600 -0.3840 L1 = 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) 运行后输出的结果与方法1 相同 .( 五) DECONV 函数调用格式 ：Q,R =dec

3、onv (B,A)( 六) roots(poly(1:n)命令调用格式 ：roots(poly(1:n)( 七) det(a*eye(size (A) - A)命令调用格式 ：b=det(a*eye(size (A) - A) 6.2 拉格朗日（ Lagrange）插值及其 MATLAB 程序6.2.1 线性插值及其 MATLAB 程序例 6.2.1 已知函数)(xf在3 , 1上具有二阶连续导数，5)( xf，且满足条件2)3(, 1)1 (ff.求线性插值多项式和函数值)5.1 (f，并估计其误差. 解输入程序 X=1,3;Y=1,2; l01= poly(X(2)/( X(1)- X(2

4、), 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=1.5; Y = polyval(P,x) 运行后输出基函数l0和 l1及其插值多项式的系数向量P（略） 、插值多项式L 和插值 Y为第六章函数的插值方法精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 1 页，共 29 页学习必备欢迎下载l0 = l1 = L = Y = -1/2*x+3/2 1/2*x-1/2 1/2*x+1

5、/2 1.2500 输入程序 M=5;R1=M*abs(x-X(1)* (x-X(2)/2 运行后输出误差限为 R1 = 1.8750 例 6.2.2求函数)(xfex在 1 ,0上线性插值多项式，并估计其误差. 解输入程序 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), 运行后输出基函数l0和 l1及其插值多项式的系数向量P 和

6、插值多项式L 为l0 = l1 = P = -x+1 x -0.6321 1.0000 L = -1423408956596761/2251799813685248*x+1 输入程序 M=1;x=0:0.001:1; R1=M*max(abs(x-X(1).*(x-X(2)./2 运行后输出误差限为 R1 = 0.1250. 6.2.2 抛物线插值及其MATLAB 程序例6.2.3求将区间0, /2 分成n等份)2, 1(n，用xxfycos)(产生1n个节点，然后根据（6.9）和（ 6.13）式分别作线性插值函数)(1xP和抛物线插值函数)(2xP.用它们分别计算cos ( /6) (取四位

7、有效数字)，并估计其误差. 解输入程序 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 (l11), P = l01* Y(1)+ l11* Y(2), L=poly2sym (P),x=pi/6; Y = polyval(P,x) 运行后输出基函数l0和 l1及其插值多项式的系数向量P、插值多项式和插值为l0 = -5734161139222659/9007199254740992*x+1 l1 = 5734161139

8、222659/9007199254740992*x P = -0.6366 1.0000 L = -5734161139222659/9007199254740992*x+1 Y = 0.6667 输入程序 M=1;x=pi/6; R1=M*abs(x-X(1)*(x-X(2)/2 运行后输出误差限为R1 = 0.2742. （2） 输入程序 X=0:pi/4:pi/2; Y =cos(X) , 精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 2 页，共 29 页学习必备欢迎下载l01= conv (poly(X(2), poly(X(3)/( X

9、(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), P = l01* Y(1)+ l11* Y(2) + l21* Y(3), L=poly2sym (P),x=pi/6; Y = polyval(P,x) 运行后输出基函数l01、l11和 l21及其

10、插值多项式的系数向量P、插值多项式 L和插值 Y为l0 = 228155022448185/281474976710656*x2-2150310427208497/1125899906842624*x+1 l1 = -228155022448185/140737488355328*x2+5734161139222659/2251799813685248*x l2 = 228155022448185/281474976710656*x2-5734161139222659/9007199254740992*x P = -0.3357 -0.1092 1.0000 L= -6048313895780

11、875/18014398509481984*x2-7870612110600739/72057594037927936*x+1 Y = 0.8508 输入程序 M=1;x=pi/6; R2=M*abs(x-X(1)*(x-X(2) *(x-X(3)/6 运行后输出误差限为R2 = 0.0239. 6.2.3 n次拉格朗日（ Lagrange）插值及其 MATLAB 程序例6.2.4给出节点数据00.17)00.2(f，00.1)00.0(f，00.2)00.1(f，00.17)00.2(f，作三次拉格朗日插值多项式计算)6 .0(f，并估计其误差.解输入程序 X=-2,0,1,2; Y =17

12、,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(3)- X(2) * ( X(3)- X(4), l31= conv ( c

13、onv (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), 运行后输出基函数l0，l1，l2和l3及其插值多项式的系数向量P（略）为l0 = -1/24*x3+1/8*x2-1/12*x，l1 =1/4*x3-1/4*x2-x+1 l2 = -1/3*x3+4/3*x，l3 =1/8*x

14、3+1/8*x2-1/4*x 输入程序 L=poly2sym (P),x=0.6; Y = polyval(P,x) 运行后输出插值多项式和插值为L = Y = x3+4*x2-4*x+1 0.2560. 精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 3 页，共 29 页学习必备欢迎下载输入程序 syms M; x=0.6; R3=M*abs(x-X(1)*(x-X(2) *(x-X(3) *(x-X(4)/24 运行后输出误差限为R3 = 91/2500*M 即R3)(04.0)4(f，)00.2,00.2(. 6.2.5 拉格朗日多项式和基函

15、数的MATLAB 程序求拉格朗日插值多项式和基函数的MATLAB 主程序function C, 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); end end L1(k,:)=V; l(k,:)=poly2sym (V) end C=Y*L1;L=Y*l例6.2.5 给出节点数据03.17)15.2(f，24.7)00.1(f，05.1)01.0(f，03. 2)02.1(f，06.17)03.2(f，05.23)25. 3

16、(f， 作五次拉格朗日插值多项式和基函数，并写出估计其误差的公式.解在 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; C, L ,L1,l= lagran1(X,Y) 运行后输出五次拉格朗日插值多项式L 及其系数向量C，基函数 l 及其系数矩阵L1如下C = -0.2169 0.0648 2.1076 3.3960 -4.5745 1.0954 L = 1.0954-4.5745*x+3.3960*x2+2.1076*x3+0.0648*x4-0.2169*x5 L1

17、 = -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.3166 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.0002 l = -0.0056*x5+0.0299*x4-0.

18、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-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

19、+0.0220*x-0.0002 估计其误差的公式为)(5xR)25. 3)(03. 2)(02. 1)(01.0()00. 1)(15. 2(! 6)()6(xxxxxxf，)3.25,-2.15(.精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 4 页，共 29 页学习必备欢迎下载6.2.6 拉格朗日插值及其误差估计的MATLAB 程序拉格朗日插值及其误差估计的MATLAB 主程序function y,R=lagranzi(X,Y,x,M) n=length(X); m=length(x); for i=1:m z=x(i);s=0.0; fo

20、r k=1:n p=1.0; q1=1.0; c1=1.0; for j=1:n if j=k p=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; end R=M*q1/c1; 例 6.2.6已知5 .030sin，1707.045sin，0866.060sin，用拉格朗日插值及其误差估计的MATLAB 主程序求40sin的近似值，并估计其误差.解在MATLAB 工作窗口输入程序 x=2*pi/9; M=1; X=pi/6 ,pi/4, pi/3; Y=0.5,0.7071,0.8

21、660; y,R=lagranzi(X,Y,x,M) 运行后输出插值y 及其误差限R 为y = R = 0.6434 8.8610e-004. 6.3 牛顿（ Newton）插值及其 MATLAB 程序6.3.3 牛顿插值多项式、差商和误差公式的MATLAB 程序求牛顿插值多项式和差商的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

22、,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:1 C=conv(C,poly(X(k); d=length(C); C(d)=C(d)+A(k,k); end L(k,:)=poly2sym(C); Q=poly2sym(q1); syms M wcgs=M*Q/c1; Cw=q1/c1; 例 6.3.3 给 出 节 点 数 据03.17)15.2(f，24.7)00.1(f，05.1)01

23、.0(f，03. 2)02.1(f，06.17)03.2(f，05.23)25. 3(f,作五阶牛顿插值多项式和差商，并写出其估计误差的公式.解 （1）保存名为newpoly.m的 M 文件 . （2）输入 MATLAB 程序精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 5 页，共 29 页学习必备欢迎下载 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) 运行后输出差商矩阵A,五阶牛顿插值多项式L

24、及其系数向量C, 插值余项公式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 23.0500 4.9098 -4.4715 -3.5056 -1.0867 -0.2169 C = -0.2169 0.0648 2.1076 3.3960 -4.5745 1.0954 L = -7813219284746629/3602879701896

25、3968*x5+583849564517807/9007199254740992*x4+593245028711263/281474976710656*x3+3823593773002357/1125899906842624*x2-321902673270315/70368744177664*x+308328649211299/281474976710656 wcgs = 1/720*M*(x6-79/25*x5-14201/2500*x4+4934097026900981/281474976710656*x3+154500712237335/35184372088832*x2-8170642

26、380559269/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* x5.估计其误差的公式为)(5xR)25. 3)(03. 2)(02. 1)(01. 0()00. 1)(15. 2(! 6)()6(xxxxxxf，)3.25,-2.15(. 例 6.3.4 求函数7)(xfe5/x在6,2上五阶牛顿

27、插值多项式，估计其误差的公式和误差限公式.用它们计算)156.3(f，并估计其误差. 解 （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), 运行后输出差商矩阵A, 五阶牛顿插值多项式L 及其系数向量C, 插值余项公式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 0 0 -2.9035 0.6297 -0.06

28、83 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.0000 C = 0.0000 -0.0004 0.0089 -0.1389 1.3985 -6.9991 L = 9721799720875/1152921504606846976*x5-3503994098647815/9223372036854775808*x4+160742008798419/18014398509481984*x3-1251152213853501/9007199254740992*

29、x2+6298131904328647/4503599627370496*x-3940156929554013/562949953421312 wcgs = 1/720*M*(x6-24*x5+1172/5*x4-5952/5*x3+7276634802928539/2199023255552*x2-5237461186650519/1099511627776*x+6085939356447121/219902精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 6 页，共 29 页学习必备欢迎下载3255552) Cw = 0.0014 -0.0333

30、 0.3256 -1.6533 4.5959 -6.6159 3.8438 M = -1.3494e-004 （2）输入 MATLAB 程序 syms x wcgs1=1/720*M*1/720*M*(x6-24*x5+1172/5*x4-5952/5*x3+7276634802928539/2199023255552*x2-5237461186650519/1099511627776*x+6085939356447121/2199023255552), 运行后输出误差限公式wcgs1 如下wcgs1 = 5565367633581443/158456325028528675187087900

31、672*x6-16696102900744329/19807040628566084398385987584*x5+1630652716639362799/198070406285660843983859875840*x4-517579189923074199/12379400392853802748991242240*x3+40497147813610772932297365501777/348449143727040986586495598010130648530944*x2-29148396970323855270001164718917/174224571863520493293247

32、799005065324265472*x+33870489914310283926665276375603/348449143727040986586495598010130648530944 （3）输入 MATLAB 程序 x=3.156; y=9721799720875/1152921504606846976*x5-3503994098647815/9223372036854775808*x4+160742008798419/18014398509481984*x3-1251152213853501/9007199254740992*x2+6298131904328647/45035996

33、27370496*x-3940156929554013/562949953421312 wcgs2=1/720*M*(x6-24*x5+1172/5*x4-5952/5*x3+7276634802928539/2199023255552*x2-5237461186650519/1099511627776*x+6085939356447121/2199023255552) 运行后输出)156.3(f的近似值 y，及其误差限wcgs 2 如下y = -3.7237 wcgs2 = -2.4764e-007 6.3.4 牛顿插值及其误差估计的MATLAB 程序牛顿插值及其误差估计的MATLAB 主程

34、序function y,R= newcz(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)=(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:1 C=conv(C,poly(X(k);d=lengt

35、h(C); C(d)=C(d)+A(k,k); end y(k)= polyval(C, z); end R=M*q1/c1; 例6.3.5 已知5 .030sin，1707.045sin，0866.060sin，用牛顿插值精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 7 页，共 29 页学习必备欢迎下载法求40sin的近似值，估计其误差，并与例 6.2.6的计算结果比较.解方法 1（牛顿插值及其误差估计的MATLAB 主程序） 输入 MATLAB 程序 x=2*pi/9;M=1; X=pi/6 ,pi/4, pi/3; Y=0.5,0.7071

36、,0.8660; y,R= newcz(X,Y,x,M) 运行后输出y = R = 0.6434 8.8610e-004 方法 2 （求牛顿插值多项式和差商的MATLAB 主程序） 输入 MATLAB 程序 x=2*pi/9; X=pi/6 ,pi/4, pi/3; Y=0.5,0.7071,0.8660; M=1; A,C,L,wcgs,Cw= newploy(X,Y), y=polyval(C,x) 运行后输出结果A = 0.5000 0 0 0.7071 0.7911 0 0.8660 0.6070 -0.3516 C = -0.3516 1.2513 -0.0588 L = -1583

37、578379808357/4503599627370496*x2+1408883391907005/1125899906842624*x-132405829044691/2251799813685248 wcgs = 1/6*M*(x3-3/4*x2*pi+4012734077357799/2251799813685248*x-7757769783530263/18014398509481984) Cw = 0.1667 -0.3927 0.2970 -0.0718 y = 0.6434 上述两种方法计算y 的结果相同 .6.3.5 牛顿插值法的 MATLAB 综合程序求牛顿插值多项式、差商、

38、插值及其误差估计的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)=(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:1 C

39、=conv(C,poly(X(k); d=length(C);C(d)=C(d)+A(k,k); end y(k)= polyval(C, z); end R=M*q1/c1;L(k,:)=poly2sym(C); 例6.3.6给出节点数据00.27)00.4(f，00.1)00.0(f，00.2)00.1(f，00.17)00.2(f，作三阶牛顿插值多项式，计算)345.2(f，并估计其误差. 解 首先将名为 newdscg.m 的程序保存为 M文件， 然后在 MA TLAB 工作窗口输入程序 syms M,X=-4,0,1,2; Y =27,1,2,17; x=-2.345; y,R,A,

40、C,P=newdscg(X,Y,x,M) 精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 8 页，共 29 页学习必备欢迎下载运行后输出插值y)345.2(f及其误差限公式R，三阶牛顿插值多项式P 及其系数向量C，差商的矩阵A 如下y = 22.3211 R = 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.9167 C = 0.9

41、167 4.2500 -4.1667 1.0000 P = 11/12*x3+17/4*x2-25/6*x+1 例 6.3.7 将区间0, /2 分成n等份)3, 2(n，用xxfysin)(产生1n个节点，求二阶和三阶牛顿插值多项式)(2xP和)(3xP.用它们分别计算sin (/7) (取四位有效数字 )，并估计其误差. 解首先将名为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=newds

42、cg(X1,Y1,x,M), y2,R2,A2,C2,P3=newdscg(X2,Y2,x,M) 运行后输出插值y1=)7/sin()7/(2P和 y2=)7/sin()7/(3P及其误差限 R1和 R2，二阶和三阶牛顿插值多项式P2和 P3及其系数向量C1和 C2，差商的矩阵A1和 A2如下y1 = 0.4548 R1 = 0.0282 A1 = 0 0 0 0.7071 0.9003 0 1.0000 0.3729 -0.3357 C1 = -0.3357 1.1640 0 P2 = -3024156947890437/9007199254740992*x2+163820246322191

43、/140737488355328*x y2 = 0.4345 R2 = 9.3913e-004 A2 = 0 0 0 0 0.5000 0.9549 0 0 0.8660 0.6991 -0.2443 0 1.0000 0.2559 -0.4232 -0.1139 C2 = -0.1139 -0.0655 1.0204 0 P3 = -1025666884451963/9007199254740992*x3-4717668559161127/72057594037927936*x2+4595602396951547/4503599627370496*x6.4 埃尔米特（ Hermite ）插值

44、及其 MATLAB 程序精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 9 页，共 29 页学习必备欢迎下载6.4.3 埃尔米特插值多项式和误差公式的MATLAB 程序求埃尔米特插值多项式和误差公式的MATLAB 主程序function Hc, Hk,wcgs,Cw= hermite (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+1 if k=i s=s+(1/(X(k)-X(i); V=conv(

45、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); end for i=1:t c1=c1*i; end syms M,wcgs=M*Q/c1; Cw=q/c1; 例 6.4.3 给定函数)(xf在点2/, 4/, 6/210xxx处的函数值5 . 0)(0xf, 1707. 0)(

46、1xf，1)(2xf和导数值0866.0)(0xf，1707.0)(1xf，0000.0)(0xf，且1)()6(xf，求函数)(xf在点210,xxx处的五阶埃尔米特插值多项式)(5xH和误差公式，计算)1.567(f并估计其误差. 解 （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) 运行后输出五阶埃尔米特插值多项式Hk及其系数向量Hc，误差公式 wcgs及其系数向量 Cw如下Hc

47、 = 1.0e+003 * 0.1912 -0.9273 1.6903 -1.4380 0.5751 -0.0866 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+7446

48、708432019761/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 当Mxf1)()6(时的误差公式为R=0.001 4* x6-0.008 0* x5+0.018 4* x4-0.021 5* x3+0.013

49、 6* x2-0.004 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-6094057精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 10 页，共 29 页学习必

50、备欢迎下载839958843/70368744177664, wcgs=1/720*M*(x6-11/6*x5*pi+7446708432019761/562949953421312*x4-4363745503235773/281474976710656*x3+21569239021155/2199023255552*x2-7178073637328281/2251799813685248*x+3758430567659515/9007199254740992) 运行后输出)1.567(f的近似值Hk 及其误差wcgs 如下Hk = 2.5265 wcgs = 1.3313e-008 6.5

51、分段插值及其 MATLAB 程序6.5.1 高次插值的振荡和MATLAB 程序例 6.5.1作 下 列 函 数 在 指 定 区 间 上 的n次 拉 格 朗 日 插 值 多 项 式)(xLn)10,8,6,4,2(n的图形，并讨论插值多项式)(xLn的次数与误差)(xRn的关系 . （1）)432sin3tan(cos()(2xxxf，,x；（2）211)(xxf，5,5x. 解将计算n次拉格朗日插值多项式)(xLn的值的 MATLAB 程序保存名为lagr1.m 的M文件. function y=lagr1(x0,y0,x) n=length(x0); m=length(x); for i=1

52、:m z=x(i);s=0.0; for k=1:n p=1.0; for j=1:n if j=k p=p*(z-x0(j)/(x0(k)-x0(j); end end s=p*y0(k)+s; end y(i)=s; end （1）现提供作n次拉格朗日插值多项式)(xLn的图形的 MATLAB 程序，保存名为Li1.m 的M文件m=150; x=-pi:2*pi/(m-1): pi; y=tan(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),pause n=3; x

53、0=-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 Li1.m 回车运行后，便会逐次画出)(xf在区间,上的n次拉格朗日插值多项式)(xLn)10,8,6,4,2(n的图形 (略). （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),pause n=3; x0=-5:10/(n-1):5; y0=1

54、./(1+x0.2); y1=lagr1(x0,y0,x);hold on, plot(x,y1,g),gtext(n=2),pause,hold off n=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 off n=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 off n=9; x0=-5:

55、10/(n-1):5; y0=1./(1+x0.2); y4=lagr1(x0,y0,x);hold on, plot(x,y4,r:),gtext(n=8),pause,hold off n=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( 高次拉格朗日插值的振荡 ) 回车运行后，便会逐次画出)1/(12xy在区间5,5上的n次拉格朗日插值多项式)(xLn)10,8,6,4,2(n的图形 (略).6.5.4 作有关分段线性插值图形的MATLAB 程

56、序作有关分段线性插值图形的MATLAB主程序function s= xxczhjt1(x0,y0,xi,x,y) s= interp1(x0,y0,xi); Sn= interp1(x0,y0,x0); plot(x0,y0,o,x0,Sn,-,xi,s,*,x,y,-.) legend( 节点 (xi,yi), 分段线性插值函数Sn(x), 插值点 (x,s), 被插值函数 y ) 我们也可以直接在在MATLAB 工作窗口编程序. 例如， x0 =-6:6; y0 =sin(x0); xi = -6:.25:6; yi = interp1(x0,y0,xi); x=-6:0.001:6; y

57、=sin(x);plot(x0,y0,o,xi,yi,x,y,:), legend( 节点 (xi,yi), 分段线性插值函数 , 被插值函数 y=sinx ) title(y=sinx及其分段线性插值函数和节点的图形) x0 =-6:6; y0 =cos(x0); xi = -6:.25:6; yi = interp1(x0,y0,xi); 精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 12 页，共 29 页学习必备欢迎下载x=-6:0.001:6; y=cos(x); plot(x0,y0,o,xi,yi,x,y,:), legend( 节点

58、 (xi,yi), 分段线性插值函数 , 被插值函数 y=cosx ) title(y=cosx及其分段线性插值函数和节点的图形)例6.5.5设函数22511)(xxf， 在区间 1 , 1上取等距节点),(iiyx10, 2, 1, 0i, 构造分段线性插值函数)(xSn，用 MATLAB 程序计算各小区间的中点ix处)(xSn的值，作出节点，插值点，)(xf和)(xSn的图形 . 解节点的横坐标和插值点等取值与例6.5.4相同 . 在 MATLAB 工作窗口输入程序x0=-1:0.2:1; y0=1./(1+25.*x0.2); xi=-0.9:0.2:0.9; b=max(x0); a=

59、min(x0);x=a:0.001:b; y=1./(1+25.*x.2); s=xxczhjt1(x0,y0,xi,x,y), title(y=1/(1+25 x2)的分段线性插值的有关图形) 运行后屏幕显示各小区间中点ix处)(xSn的值，出现节点，插值点，)(xf和)(xSn的图形(略)s = Columns 1 through 4 0.04864253393665 0.07941176470588 0.15000000000000 0.35000000000000 Columns 5 through 8 0.75000000000000 0.75000000000000 0.35000

60、000000000 0.15000000000000 Columns 9 through 10 0.07941176470588 0.04864253393665 6.5.5 用 MATLAB 计算有关分段线性插值的误差例6.5.8设函数)432sin3tan(cos()(2xxxf,在区间,上取等距节点),(iiyx，9, 2, 1, 0i，构造分段线性插值函数)(xSn. (1)用MATLAB 程序计算各小区间中点ix处)(xSn的值，作出节点， 插值点，)(xf和)(xSn的图形；(2) 用MATLAB 程序计算各小区间中点处)(xSn的值及其相对误差；(3) 用MATLAB 程序估计)

61、(maxxfx和)(xSn在区间,上的误差限 . 解在 MATLAB 工作窗口输入程序h=2*pi/9; x0=-pi:h:pi; y0=tan(cos(sqrt(3)+sin(2*x0)./(3+4*x0.2); xi=-pi+h/2:h:pi-h/2; fi=tan(cos(3(1/2)+sin(2*xi)./(3+4*xi.2); b=max(x0); a=min(x0);x=a:0.001:b; y=tan(cos(sqrt(3)+sin(2.*x)./(3+4*x.2); si=xxczhjt1(x0,y0,xi,x,y); Ri=abs(fi-si)./fi); xi,fi,si,

62、Ri,i=xi,fi,si,Ri title(y=tan(cos(sqrt(3)+sin(2 x)/(3+4 x2)的分段线性插值的有关图形 ) 运行后屏幕显示Ri（略），并且作出节点，插值点，)(xf和)(xSn的图形 (略). 在 MATLAB 工作窗口输入程序精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 13 页，共 29 页学习必备欢迎下载 syms x y=tan(cos(3(1/2)+sin(2*x)/(3+4*x2);yxx=diff(y,x,2), 运行后屏幕显示函数)(xf的二阶导数)( xf（略） . 在 MATLAB 工作窗

63、口输入程序 h=2*pi/9; x=-pi:0.0001:-pi; yxx=2*tan(cos(3(1/2)+sin(2.*x)./(3+4.*x.2).*(1+tan(cos(3.(1/2)+sin(2.*x)./(3+4.*x.2).2).*sin(3(1/2)+sin(2.*x)./(3+4.*x.2).2.*(2*cos(2.*x)./(3+4*x.2)-8*(3(1/2)+sin(2.*x)./(3+4.*x.2).2.*x).2-(1+tan(cos(3(1/2)+sin(2.*x)./(3+4.*x.2).2).*cos(3(1/2)+sin(2.*x)./(3+4.*x.2).

64、*(2.*cos(2.*x)./(3+4.*x.2)-8*(3(1/2)+sin(2.*x)./(3+4.*x.2).2.*x).2-(1+tan(cos(3(1/2)+sin(2.*x)./(3+4.*x.2).2).*sin(3.(1/2)+sin(2.*x)./(3+4.*x.2).*(-4.*sin(2.*x)./(3+4.*x.2)-32*cos(2.*x)./(3+4.*x.2).2.*x+128*(3(1/2)+sin(2.*x)./(3+4.*x.2).3.*x.2-8*(3(1/2)+sin(2.*x)./(3+4.*x.2).2); myxx=max(yxx), R=(h2

65、)* myxx/8 运行后屏幕显示myxx =)(maxxfx和)(xSn在区间,上的误差限R如下myxx = R = -0.02788637150664 -0.00169893490711 6.6 分段埃尔米特插值及其MATLAB 程序6.6.2 分段埃尔米特插值的MATLAB 程序调用格式一：YI=interp1(X,Y,XI,pchip) 调用格式二 ：YI=interp1(X,XI,pchip)例 6.6.5试用 MA TLAB 程序计算例6.6.1中在各小区间中点处分段三次埃尔米特插值)(2/1inxH及其相对误差 . 解在 MATLAB 工作窗口输入程序 h=0.2;x0=-1:h

66、:1;y0=1./(1+25.*x0.2); xi=-0.9:h:0.9; fi=1./(1+25.*xi.2); yi=interp1(x0,y0,xi,pchip); Ri=abs(fi-yi)./fi); xi,fi,yi,Ri,i=xi,fi,yi,Ri 运行后屏幕显示各小区间中点xi处的函数值fi，插值 si，相对误差值Ri (略). 6.6.3 作有关分段埃尔米特插值图形的MATLAB 程序作有关分段埃尔米特插值图形的MATLAB 主程序function H=hermitetx(x0,y0,xi,x,y) H= interp1(x0,y0,xi,pchip); Hn= interp

67、1(x0,y0,x,pchip); plot(x0,y0,o,x,Hn,-,xi,H,*,x,y,-.) legend( 节点 (xi,yi), 分段埃尔米特插值函数 , 插值点 (x,H), 被插值函数 y ) 我们也可以直接在在MATLAB工作窗口编程序，例如， x0 =-6:6; y0 =sin(x0); xi = -6:.25:6; yi = interp1(x0,y0,xi,pchip); x=-6:0.001:6; y=sin(x); plot(x0,y0,o,xi,yi,x,y,:), legend( 节点 (xi,yi), 分段埃尔米特插值函数 , 被插值函数y=sinx) t

68、itle( y=sinx及其分段埃尔米特插值函数和节点的图形 ) x0 =-6:6; y0 =cos(x0); xi = -6:.25:6;yi = interp1(x0,y0,xi,pchip); x=-6:0.001:6; y=cos(x); plot(x0,y0,o,xi,yi,x,y,:), legend( 节点 (xi,yi), 分段埃尔米特插值函数 , 被插值函数y=cosx) 精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 14 页，共 29 页学习必备欢迎下载title( y=cosx及其分段埃尔米特插值函数和节点的图形) 例6.6

69、.6设函数211)(xxf定义在区间5, 5上，节点 ( X(i),f(X (i)的横坐标向量 X的元素是首项a=-5，末项 b=5，公差 h=1.5的等差数列，构造三次分段埃尔米特插值函数)(3,xHn.把区间5,5分成 20等份，构成 20个小区间， 用MATLAB 程序计算各小区间中点ix处)(3,xHn的值，并作出节点，插值点，)(xf和)(3,xHn的图形 .解在 MATLAB 工作窗口输入程序x0=-5:1.5:5; y0=1./(1+x0.2); x1=-4.75:0.5:4.75; x=-5:0.001:5; y=1./(1+x.2); H= hermitetx(x0,y0,x

70、1,x,y) title( 函数 y=1/(1+x2)及其分段埃尔米特插值函数，插值，节点 (xi,yi) 的图形 )运行后屏幕显示各小区间中点ix处)(3,xHn的值，出现节点，插值点，)(xf和)(3,xHn的图形 (略).例6.6.7设函数xxxfcos5.0)(定义在区间,上，取7n，按等距节点构造分段埃尔米特插值函数)(3,7xH，用 MATLAB 程序计算各小区间中点ix处)(3,7xH的值，作出节点，插值点，)(xf和)(3, 7xH的图形 .解记节点的横 坐标7,2 ,1 ,0,7/2,ihihxi插 值点)(21121iiixxx，6,2, 1 ,0i.在 MATLAB工作窗

71、口输入程序 h=2*pi/7; x0=-pi:h:pi; y0=0.5.*x0-cos(x0); xi=-pi+h/2:h:pi-h/2; b=max(x0); a=min(x0); x=a:0.001:b; y=0.5.*x-cos(x); H= hermitetx(x0,y0,xi,x,y) title( 函数y=0.5x-cos(x)及其分段埃尔米特插值函数，插值，节点(xi,yi) 的图形 )运行后屏幕显示各小区间中点ix处)(3,7xH的值，出现节点， 插值点，)(xf和)(3,7xH的图形 (略).6.6.4 用 MATLAB 计算有关分段埃尔米特插值的误差例6.6.8 设函数22

72、511)(xxf定义在区间 1,1上，取10n，按等距节点构造分段埃尔米特插值函数)(3,xHn，用 MATLAB 程序在 1,1上计算)(max)4(11xfx和)(3,xHn的误差公式和误差限.解在 MATLAB 工作窗口输入程序 syms h,x=-1:0.0001:1; yxxxx=150000000./(1+25.*x.2).5.*x.4-4500000./(1+25.*x.2).4.*x.2+15000./(1+25.*x.2).3; myxxxx=max(yxxxx), R=(h4)* abs(myxxxx/384) 运行后输出)(xf的 4 阶导数在区间 1, 1上绝对值的最大

73、值myxxxx 和)(3,xHn在区间 1, 1上的误差公式myxxxx 为myxxxx = R = 15000 625/16*h4 （4）在 MATLAB 工作窗口输入程序精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 15 页，共 29 页学习必备欢迎下载 h=0.2; R =625/16*h4 运行后输出误差限为R = 0.06250000000000 例6.6.9 设函数)432sin3tan(cos()(2xxxf定义在区间,上，取9n，按等距节点构造分段埃尔米特插值函数)(3,xHn. (1)用MATLAB 程序计算各小区间中点ix处)

74、(3,xHn的值， 作出节点， 插值点，)(xf和)(3,xHn的图形；(2) 并用 MATLAB 程序计算各小区间中点处)(3,xHn的值及其相对误差；(3) 用MATLAB 程序求)(max)4(xfx和)(3,xHn在区间,上的误差公式和各插值的误差限.解 （1）记节点的横坐标9 , 2, 1 , 0, 9/2,ihihxi，插值点)(21121iiixxx，8 , 2 , 1 , 0i.在 MATLAB 工作窗口输入程序h=2*pi/9; x0=-pi:h:pi; y0=tan(cos(sqrt(3)+sin(2*x0)./(3+4*x0.2); xi=-pi+h/2:h:pi-h/2

75、; fi=tan(cos(3(1/2)+sin(2*xi)./(3+4*xi.2); b=max(x0); a=min(x0); x=a:0.001:b; y=tan(cos(3(1/2)+sin(2.*x)./(3+4*x.2); Hi= hermitetx(x0,y0,xi,x,y); Ri=abs(fi-Hi)./fi); xi,fi,Hi,Ri,i=xi,fi,Hi,Ri title( 函数 y= tan(cos(sqrt(3)+sin(2x)/(3+4x2)及其分段埃尔米特插值函数，插值，节点(xi,yi) 的图形 )运行后屏幕显示各小区间中点xi处的函数值fi，插值 Hi，相对误差

76、值Ri，并且作出节点，插值点，)(xf和)(3,xHn的图形 (略). （2）在 MATLAB 工作窗口输入程序 syms x y=tan(cos(3(1/2)+sin(2*x)/(3+4*x2); yxxxx=diff(y,x,4),%simple(yxxxx) 运 行 后 屏 幕 显 示 函 数)(xf的4 阶 导 数)()4(xf， 然 后 将 输 出 的)()4(xf编 程 求)(max)4(xfx和)(3,xHn及其在区间,上的误差限的MA TLAB 程序如下syms h,x=-pi:0.0001:pi; yxxxx=-12.*(1.+tan(cos(3.(1./2)+sin(2.*

77、x)./(3.+4.*x.2).2).2.*sin(3.(1./2)+sin(2.*x)./(3.+4.*x.2).3.*(2.*cos(2.*x)./(3.+4.*x.2)-8.*(3.(1./2)+sin(2.*x)./(3.+4.*x.2).2.*x).2.*(-4.*sin(2.*x)./(3.+4.*x.2)-32.*cos(2.*x)./(3.+4.*x.2).2.*x+128.*(3.(1./2)+sin(2.*x)./(3.+4.*x.2).3.*x.2.-8.*(3.(1./2)+sin(2.*x)./(3.+4.*x.2).2)+16.*(1.+tan(cos(3.(1./

78、2)+sin(2.*x)./(3.+4.*x.2).2).2.*sin(3.(1./2)+sin(2.*x)./(3.+4.*x.2).4.*(2.*cos(2.*x)./(3.+4.*x.2)-8.*(3.(1./2)+sin(2.*x)./(3.+4.*x.2).2.*x).4.*tan(cos(3.(1./2)+sin(2.*x)./(3.+4.*x.2)+8.*tan(cos(3.(1./2)+sin(2.*x)./(3.+4.*x.2).3.*(1.+tan(cos(3.(1./2)+sin(2.*x)./(3.+4.*x.2).2).*sin(3.(1./2)+sin(2.*x).

79、/(3+4.*x.2).4.*(2.*cos(2.*x)./(3+4.*x.2)-8.*(3.(1./2)+sin(2.*x)./(3.+4.*x.2).2.*x).4-8.*tan(cos(3.(1./2)+sin(2.*x)./(3.+4.*x.2).*(1.+tan(cos(3.(1./2)+sin(2.*x)./(3.+4*x.2).2).*sin(3.(1./2)+sin(2.*x)./(3.+4.*x.2).2.*(2.*cos(2.*x)./(3.+4.*x.2)-8*(3.(1./2)+sin精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - -

80、 -第 16 页，共 29 页学习必备欢迎下载(2.*x)./(3.+4*x.2).2.*x).4+6.*(1+tan(cos(3.(1./2)+sin(2.*x)./(3+4.*x.2).2).*sin(3.(1./2)+sin(2.*x)./(3+4.*x.2).*(2.*cos(2.*x)./(3+4.*x.2)-8.*(3.(1./2)+sin(2.*x)./(3+4.*x.2).2.*x).2.*(-4.*sin(2.*x)./(3+4.*x.2)-32.*cos(2.*x)./(3+4.*x.2).2.*x+128.*(3.(1./2)+sin(2.*x)./(3+4.*x.2).

81、3.*x.2-8.*(3.(1./2)+sin(2.*x)./(3+4.*x.2).2)+(1+tan(cos(3.(1./2)+sin(2.*x)./(3+4.*x.2).2).*cos(3.(1./2)+sin(2.*x)./(3+4.*x.2).*(2.*cos(2.*x)./(3+4.*x.2)-8.*(3.(1./2)+sin(2.*x)./(3+4.*x.2).2.*x).4-3.*(1+tan(cos(3.(1./2)+sin(2.*x)./(3+4.*x.2).2).*cos(3.(1./2)+sin(2.*x)./(3+4.*x.2).*(-4.*sin(2.*x)./(3+

82、4.*x.2)-32.*cos(2.*x)./(3+4.*x.2).2.*x+128.*(3.(1./2)+sin(2.*x)./(3+4.*x.2).3.*x.2-8.*(3.(1./2)+sin(2.*x)./(3+4.*x.2).2).2-4.*(1+tan(cos(3.(1./2)+sin(2.*x)./(3+4.*x.2).2).*cos(3.(1./2)+sin(2.*x)./(3+4.*x.2).*(2.*cos(2.*x)./(3+4.*x.2)-8.*(3.(1./2)+sin(2.*x)./(3+4.*x.2).2.*x).*(-8.*cos(2.*x)./(3+4.*x.

83、2)+96.*sin(2.*x)./(3+4.*x.2).2.*x+768.*cos(2.*x)./(3+4.*x.2).3.*x.2-48.*cos(2.*x)./(3+4.*x.2).2-3072.*(3.(1./2)+sin(2.*x)./(3+4.*x.2).4.*x.3+384.*(3.(1./2)+sin(2.*x)./(3+4.*x.2).3.*x)-(1+tan(cos(3.(1./2)+sin(2.*x)./(3+4.*x.2).2).*sin(3.(1./2)+sin(2.*x)./(3+4.*x.2).*(16.*sin(2.*x)./(3+4.*x.2)+256.*co

84、s(2.*x)./(3+4.*x.2).2.*x-3072.*sin(2.*x)./(3+4.*x.2).3.*x.2+192.*sin(2.*x)./(3+4.*x.2).2-24576.*cos(2.*x)./(3+4.*x.2).4.*x.3+3072.*cos(2.*x)./(3+4.*x.2).3.*x+98304.*(3.(1./2)+sin(2.*x)./(3+4.*x.2).5.*x.4-18432.*(3.(1./2)+sin(2.*x)./(3+4.*x.2).4.*x.2+384.*(3.(1./2)+sin(2.*x)./(3+4.*x.2).3)-12.*(1+tan

85、(cos(3.(1./2)+sin(2.*x)./(3+4.*x.2).2).2*sin(3.(1./2)+sin(2.*x)./(3+4.*x.2).2.*(2.*cos(2.*x)./(3+4.*x.2)-8.*(3.(1./2)+sin(2.*x)./(3+4.*x.2).2.*x).4.*cos(3.(1./2)+sin(2.*x)./(3+4.*x.2)-24.*tan(cos(3.(1./2)+sin(2.*x)./(3+4.*x.2).2.*(1+tan(cos(3.(1./2)+sin(2.*x)./(3+4.*x.2).2).*sin(3.(1./2)+sin(2.*x)./

86、(3+4.*x.2).3.*(2.*cos(2.*x)./(3+4.*x.2)-8.*(3.(1./2)+sin(2.*x)./(3+4.*x.2).2.*x).2.*(-4.*sin(2.*x)./(3+4.*x.2)-32.*cos(2.*x)./(3+4.*x.2).2.*x+128.*(3.(1./2)+sin(2.*x)./(3+4.*x.2).3.*x.2-8.*(3.(1./2)+sin(2.*x)./(3+4.*x.2).2)-24.*tan(cos(3.(1./2)+sin(2.*x)./(3+4.*x.2).2.*(1+tan(cos(3.(1./2)+sin(2.*x).

87、/(3+4.*x.2).2).*sin(3.(1./2)+sin(2.*x)./(3+4.*x.2).2.*(2.*cos(2.*x)./(3+4.*x.2)-8.*(3.(1./2)+sin(2.*x)./(3+4.*x.2).2.*x).4.*cos(3.(1./2)+sin(2.*x)./(3+4.*x.2)+36.*tan(cos(3.(1./2)+sin(2.*x)./(3+4.*x.2).*(1+tan(cos(3.(1./2)+sin(2.*x)./(3+4.*x.2).2).*sin(3.(1./2)+sin(2.*x)./(3+4.*x.2).*(2.*cos(2.*x)./

88、(3+4.*x.2)-8.*(3.(1./2)+sin(2.*x)./(3+4.*x.2).2.*x).2.*cos(3.(1./2)+sin(2.*x)./(3+4.*x.2).*(-4.*sin(2.*x)./(3+4.*x.2)-32.*cos(2.*x)./(3+4.*x.2).2.*x+128.*(3.(1./2)+sin(2.*x)./(3+4.*x.2).3.*x.2-8.*(3.(1./2)+sin(2.*x)./(3+4.*x.2).2)+6.*tan(cos(3.(1./2)+sin(2.*x)./(3+4.*x.2).*(1+tan(cos(3.(1./2)+sin(2.

89、*x)./(3+4.*x.2).2).*cos(3.(1./2)+sin(2.*x)./(3+4.*x.2).2.*(2.*cos(2.*x)./(3+4.*x.2)-8.*(3.(1./2)+sin(2.*x)./(3+4.*x.2).2.*x).4+6.*tan(cos(3.(1./2)+sin(2.*x)./(3+4.*x.2).*(1+tan(cos(3.(1./2)+sin(2.*x)./(3+4.*x.2).2).*sin(3.(1./2)+sin(2.*x)./(3+4.*x.2).2.*(-4.*sin(2.*x)./(3+4.*x.2)-32.*cos(2.*x)./(3+4

90、.*x.2).2.*x+128.*(3.(1./2)+sin(2.*x)./(3+4.*x.2).3.*x.2-8.*(3.(1./2)+sin(2.*x)./(3+4.*x.2).2).2+8.*tan(cos(3.(1./2)+sin(2.*x)./(3+4.*x.2).*(1+tan(cos(3.(1./2)+sin(2.*x)./(3+4.*x.2).2).*sin(3.(1./2)+sin(2.*x)./(3+4.*x.2).2.*(2.*cos(2.*x)./(3+4.*x.2)-8.*(3.(1./2)+sin(2.*x)./(3+4.*x.2).2.*x).*(-8.*cos(

91、2.*x)./(3+4.*x.2)+96.*sin(2.*x)./(3+4.*x.2).2.*x+768.*cos(2.*x)./(3+4.*精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 17 页，共 29 页学习必备欢迎下载x.2).3.*x.2-48.*cos(2.*x)./(3+4.*x.2).2-3072.*(3.(1./2)+sin(2.*x)./(3+4.*x.2).4.*x.3+384.*(3.(1./2)+sin(2.*x)./(3+4.*x.2).3.*x) myxx=max(yxxxx), R=(h.4).* abs(myxx

92、./384) 将其保存为名为myxx.m的 M文件，然后在MATLAB 工作窗口输入该文件名 myxx运 行 后 屏 幕 显 示myxx =)(max)4(xfx和)(3,xHn在 区 间,上 的 误 差 公 式)(ma x384)4(4xfhRx如下myxx = R = 73.94706841647552 1734520780029061/9007199254740992*h4 最后在 MATLAB 工作窗口输入h=2*pi/9; R =1734520780029061/9007199254740992*h4 运行后屏幕显示)(3,xHn在区间,上的误差限R = 0.045744530299

93、48 6.7 三次样条（ Spline ）及其 MATLAB 程序6.7.4 用一阶导数计算的几种样条函数和MATLAB 程序计算压紧样条的 MATLAB 主程序functionm,H,lambda,mu,D,A,dY,sk=ClampedSp (X,Y,dyo,dyn) m=length(X);A=zeros(m,m); n=m-1;H=zeros(1,n);lambda=zeros(1,n); mu=zeros(1,n);lambda(1)=0;A(1,1)=2;A(1,2)=lambda(1); D=zeros(1,n);H(1)=X(2)-X(1);mu(1)=1-lambda(1);

94、D(1)=2*dyo; for k=1:n hk=X(k+1)-X(k);H(k+1)=hk; end H=H(2:n+1); for k=1:n-1 lambdak=H(k)/(H(k)+H(k+1); lambda(k+1)=lambdak; muk=1-lambda(k+1);mu(k)=muk; dk=3*(mu(k).*(Y(k+1)-Y(k)./H(k)+(lambda(k+1).*(Y(k+2)-Y(k+1)./H(k+1); D(k+1)=dk; end D(m)=2*dyn;mu(n)=0; n;H;lambda;mu;D; for i=1:m-1 A(i,i)=2;A(m,

95、m)=2; A(i,i+1)=lambda(i); A(i+1,i)=mu(i); end dY=AD; syms x m=length(X);S=zeros(m-1,1); for k=2:m sk=Y(k-1)*(H(k-1)-2*X(k-1)+2*x)*(x-X(k)2)/(H(k-1)3)+Y(k)*(H(k-1)+2*X(k)-2*x)*(x-X(k-1)2)/(H(k-1)3)+dY(k-1)*(x-X(k-1)*(x-X(k)2)/(H(k-1)2)+dY(k)*(x-X(k)*(x-X(k-1)2)/(H(k-1)2) end例 6.7.2 用计算压紧样条的MATLAB 主程序

96、 ClampedSp.m求例 6.7.1的问题 . 解保存 ClampedSp.m为 M 文件，输入程序精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 18 页，共 29 页学习必备欢迎下载 X=0:2:6,10,Y=0,16,36,54,82,dyo=8,dyn=7, m,H,lambda,mu,D,A,dY,sk=ClampedSp(X,Y,dyo,dyn) 运行后输出压紧样条函数sk ， X 的维数m ，由1ih), 2, 1(ni、i、i和id)1, 2, 1(ni的坐标组成的向量H、lambda 和 mu、D， （6.90）式的系数矩阵A

97、，线性方程组（6.90）的解向量dY 如下sk = 2*(6-2*x)*x2+2*x*(x-2)2+9/4*(x-2)*x2 sk = 2*(-2+2*x)*(x-4)2+9/2*(10-2*x)*(x-2)2+9/4*(x-2)*(x-4)2+5/2*(x-4)*(x-2)2 sk = 9/2*(-6+2*x)*(x-6)2+27/4*(14-2*x)*(x-4)2+5/2*(x-4)*(x-6)2+2*(x-6)*(x-4)2 sk = 27/32*(-8+2*x)*(x-10)2+41/32*(24-2*x)*(x-6)2+1/2*(x-6)*(x-10)2+7/16*(x-10)*(x

98、-6)2 m = 5 H = 2 2 2 4 lambda = 0 0.5000 0.5000 0.3333 mu = 0.5000 0.5000 0.6667 0 D = 16.0000 27.0000 28.5000 25.0000 14.0000 A = 2.0000 0 0 0 0 0.5000 2.0000 0.5000 0 0 0 0.5000 2.0000 0.5000 0 0 0 0.6667 2.0000 0.3333 0 0 0 0 2.0000 dY = 8.0000 9.0000 10.0000 8.0000 7.0000输入程序 x2=3, s2=2*(-2+2*x)

99、*(x-4)2+9/2*(10-2*x)*(x-2)2+9/4*(x-2)*(x-4)2+5/2*(x-4)*(x-2)2 x4=8, s4=27/32*(-8+2*x)*(x-10)2+41/32*(24-2*x)*(x-6)2+1/2*(x-6)*(x-10)2+7/16*(x-10)*(x-6)2 运行后输出)3()3(2sf和)8()8(4sf的值如下x2 = s2 = x4 = s4 = 3 25.7500 8 68.5000 例 6.7.2和例 6.7.1计算的结果完全相同. 计算 自然样条的MATLAB 主程序function m,H,lambda,mu,D,A,dY,sk=Na

100、turalSp(X,Y) m=length(X);A=zeros(m,m); n=m-1;H=zeros(1,n);lambda=zeros(1,n); mu=zeros(1,n);lambda(1)=1;A(1,1)=2;A(1,2)=lambda(1); D=zeros(1,n);H(1)=X(2)-X(1);mu(1)=1;D(1)=3*(Y(2)-Y(1); for k=1:n hk=X(k+1)-X(k);H(k+1)=hk; 精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 19 页，共 29 页学习必备欢迎下载end H=H(2:n+1

101、); for k=1:n-1 lambdak=H(k)/(H(k)+H(k+1); lambda(k+1)=lambdak; muk=1-lambda(k+1);mu(k)=muk; dk=3*(mu(k).*(Y(k+1)-Y(k)./H(k)+(lambda(k+1).*(Y(k+2)-Y(k+1)./H(k+1); D(k+1)=dk; end D(m)= 3*(Y(m)-Y(m-1)/H(m-1);mu(n)=1; n;H;lambda;mu;D; for i=1:m-1 A(i,i)=2;A(m,m)=2; A(i,i+1)=lambda(i); A(i+1,i)=mu(i); en

102、d dY=AD; syms x m=length(X);S=zeros(m-1,1); for k=2:m sk=Y(k-1)*(H(k-1)-2*X(k-1)+2*x)*(x-X(k)2)/(H(k-1)3)+Y(k)*(H(k-1)+2*X(k)-2*x)*(x-X(k-1)2)/(H(k-1)3)+dY(k-1)*(x-X(k-1)*(x-X(k)2)/(H(k-1)2)+dY(k)*(x-X(k)*(x-X(k-1)2)/(H(k-1)2) end 例 6.7.3 用计算自然样条的MATLAB 主程序 NaturalSp.m求例 6.7.1的问题 . 解保存 ClampedSp.m为

103、M 文件，输入程序 X=0:2:6,10,Y=0,16,36,54,82,dyo=8,dyn=7, m,H,lambda,mu,D,A,dY,sk=NaturalSp(X,Y) 运行后输出自然样条函数sk ，X 的维数m ，由1ih), 2, 1(ni、i、i和id)1, 2, 1(ni的坐标组成的向量H、lambda 和 mu、D， （6.88）式的系数矩阵A，线性方程组（6.88）的解向量dY 如下sk = 2*(6-2*x)*x2+915/172*x*(x-2)2+117/86*(x-2)*x2 sk = 2*(-2+2*x)*(x-4)2+9/2*(10-2*x)*(x-2)2+117

104、/86*(x-2)*(x-4)2+471/172*(x-4)*(x-2)2 sk = 9/2*(-6+2*x)*(x-6)2+27/4*(14-2*x)*(x-4)2+471/172*(x-4)*(x-6)2+333/172*(x-6)*(x-4)2 sk = 27/32*(-8+2*x)*(x-10)2+41/32*(24-2*x)*(x-6)2+333/688*(x-6)*(x-10)2+285/688*(x-10)*(x-6)2 m = 5 H = 2 2 2 4 lambda = 1 1/2 1/2 1/3 mu = 1/2 1/2 2/3 1 D = 48 27 57/2 25 21

105、 A = 2 1 0 0 0 1/2 2 1/2 0 0 0 1/2 2 1/2 0 0 0 2/3 2 1/3 0 0 0 1 2 dY = 915/43 234/43 471/43 333/43 285/43 精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 20 页，共 29 页学习必备欢迎下载输入程序 x2=3, s2=2*(-2+2*x)*(x-4)2+9/2*(10-2*x)*(x-2)2+117/86*(x-2)*(x-4)2+471/172*(x-4)*(x-2)2 x4=8, s4=s4=27/32*(-8+2*x)*(x-10)2

106、+41/32*(24-2*x)*(x-6)2+333/688*(x-6)*(x-10)2+285/688*(x-10)*(x-6)2 运行后输出)3()3(2sf和)8()8(4sf的值如下x2 = s2 = x 4= s4 = 3 24.6221 8 68.5581 例 6.7.3和例 6.7.2用的方法不同，所以计算的结果不同，但是两种方法计算的结果相近 . 6.7.6 用 MATLAB 计算三次样条例 6.7.6 给出节点的数据如下：x-1.00 -0.54 0.13 1.12 1.89 2.06 2.54 2.82 3.50 )(xf-2.46 -5.26 -1.87 0.05 1.6

107、5 2.69 4.56 7.89 10.31 分别求在下列条件下在插值点02.01x，56.22x处的压紧三次样条插值，并显示该样条函数的有关信息：（1）端点约束条件为5) 1(nS，16.29)50. 3(nS；（2）端点约束条件为0) 1(nS，0)50.3(nS. 解 （1）输入 MATLAB 程序 X=-1.00 -0.54 0.13 1.12 1.89 2.06 2.54 2.82 3.50; Y=-2.46 -5.26 -1.87 0.05 1.65 2.69 4.56 7.89 10.31; XI=-0.02 2.56; YI= spline (X, 5,Y,29.16,XI),

108、 PP = spline (X, 5,Y,29.16) 运行后屏幕显示压紧样条分别在02.01x，56.22x处的插值和该样条函数的有关信息如下YI = -3.1058 4.7834 PP = form: pp breaks: -1 -0.5400 0.1300 1.1200 1.8900 2.0600 2.5400 2.8200 3.5000 coefs: 8x4 double pieces: 8 order: 4 dim: 1（2）因为端点约束条件为0) 1(nS，0)50.3(nS，所以输入MATLAB 程序 YI= spline (X, 0,Y,0,XI), PP= spline (X

109、, 0,Y,0) 运行后屏幕显示压紧三次样条分别在02.01x，56. 22x的插值和该样条函数的有关信息如下YI = -3.0192 4.7501 精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 21 页，共 29 页学习必备欢迎下载PP = form: pp breaks: -1 -0.5400 0.1300 1.1200 1.8900 2.0600 2.5400 2.8200 3.5000 coefs: 8x4 double pieces: 8 order: 4 dim: 1 例 6.7.7在默认端点约束条件下，用两种方法计算在下列插值点处的

110、三次样条插值. （1）给出节点的数据与例6.7.6相同 , 插值点 XI=2.56；（2） 节点 ( X(i),Y(i)的横坐标向量X 从5到5的整数，纵坐标向量Y= （ -2.36 ， 0.85 ，1.12 ，2.36 ，2.36 ，3，4，1.46 ，0.49 ， 0.06, 0） ，插值点向量XI 是从4到4的 11个等分点 . 解 （1）输入 MATLAB 程序X=-1.00 -0.54 0.13 1.12 1.89 2.06 2.54 2.82 3.50; Y=-2.46 -5.26 -1.87 0.05 1.65 2.69 4.56 7.89 10.31; XI=2.56; Y1=

111、 spline (X, Y,XI),Y2=interp1(X,Y,XI,spline) 运行后屏幕显示三次样条插值的两种结果如下Y1 = 4.730 2，Y2 =4.730 2. （2）输入 MATLAB 程序 X= -5:5; Y= -2.36 .85 1.12 2.36 2.36 3 4 1.46 .49 .06 0; XI = linspace(-4,4,11); Y1= spline (X, Y,XI), Y2=interp1(X,Y,XI,spline) 运行后屏幕显示Y1 = 0.8500 1.0203 1.8857 2.4779 2.3683 3.0000 4.0656 2.59

112、35 0.8247 0.4074 0.0600 Y2 = 0.8500 1.0203 1.8857 2.4779 2.3683 3.0000 4.0656 2.5935 0.8247 0.4074 0.0600 例 6.7.8设函数22511)(xxf定义在区间 1, 1上，取10n，按等距节点构造分段三次样条函数)(xSn.试用MATLAB程序计算在各小区间中点处分段三次样条插值)(2/1inxS及其相对误差. 解.在 MATLAB 工作窗口输入程序 h=0.2;x0=-1:h:1;y0=1./(1+25.*x0.2);xi=-0.9:h:0.9; fi=1./(1+25.*xi.2); y

113、i=spline (x0,y0,xi); Ri=abs(fi-yi)./fi); xi,fi,yi,Ri,i=xi,fi,yi,Ri 运行后屏幕显示各小区间中点xi处的函数值fi，插值 si，相对误差值Ri(略). 6.7.7 几种作三次样条有关图像的MATLAB 程序求有关分段三次样条图形的MATLAB 主程序（一）限定端点约束条件的作图程序function S=splinetx(x0,y0,x,x,y,dy1,dyn) S = spline(x0,dy1,y0,dyn,x ); Sn = spline(x0,dy1,y0,dyn,x); plot(x0,y0,o,x,Sn,-,x ,S,

114、*,x,y,-.) legend( 节点 (xi,yi), 分段三次样条函数 , 插值点 (x,S), 被插值函数 y ) 精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 22 页，共 29 页学习必备欢迎下载（二）不限定端点约束条件的作图程序function S=splinetx1(x0,y0,xi,x,y) S= interp1(x0,y0,xi, spline); Sn= interp1(x0,y0,x, spline); plot(x0,y0,o,x,Sn,-,xi,S,*,x,y,-.) legend( 节点 (xi,yi), 分段三次样

115、条函数 , 插值点 (x,S), 被插值函数 y ) （三）自由作图程序直接在 MATLAB 工作窗口编程序，例如，subplot(2,2,1),x1=-8:4/3:-4,c1=sin(x1);xx1 = -8:0.1:-4; pp1 = interp1 (x1,c1,xx1,spline ); cc1 =sin(xx1);%pp1 = spline (x1,c1,xx1); plot(x1,c1,bo,xx1,pp1,k-,xx1,cc1,r-.) subplot(2,2,2) x2=-4:4/3:-0;c2=sin(x2); xx2 = -4:0.1:0; pp2 = interp1 (x

116、2,c2,xx2,spline ); cc2=sin(xx2);plot(x2,c2,bo,xx2,pp2,k-,xx2,cc2,r-.) title(y=sinx及其三次样条插值函数, 节点 (xi,yi)的图形 ) subplot(2,1,2) x=-8:4/3:8;c=sin(x);xx = -8:0.1:8; pp = spline(x,c,xx); cc=sin(xx); plot(x,c,bo,xx,pp,k-,xx,cc,r-.) legend( 节点 (xi,yi), 三次样条插值函数 , y=sinx 的函数 ) 或 subplot(2,2,1),x1=-8:4/3:-4,c

117、1=cos(x1);xx1 = -8:0.1:-4; Y1=interp1(x1,c1,xx1, pchip); pp1 = interp1 (x1,c1,xx1,spline ); cc1 =cos(xx1);plot(x1,c1,bo,xx1,pp1,k-,xx1,Y1,r:,xx1,cc1,g-.) subplot(2,2,2) x2=-4:4/3:-0;c2= cos(x2); xx2 = -4:0.1:0; pp2 = interp1 (x2,c2,xx2,spline ); Y2=interp1(x2,c2,xx2, pchip);cc2=cos(xx2);plot(x2,c2,b

118、o,xx2,pp2,k-,xx2,Y2,r:,xx2,cc2,g-.) title(y=cosx及其三次样条插值函数, 分段三次埃尔米特插值, 节点(xi,yi)的图形 ) subplot(2,1,2) x=-8:4/3:8;c= cos(x);xx = -8:0.1:8; pp = spline(x,c,xx);Y=interp1(x,c,xx, pchip); cc= cos(xx); plot(x,c,bo,xx,pp,k-,xx,Y,r:,xx,cc,g-.) legend( 节点 (xi,yi), 三次样条插值函数 , 分段三次埃尔米特插值 , y=cosx 的函数 ) 或 n=7,

119、h=2*pi/n; x=(-2*pi+h)/2:h:(2*pi-h)/2; y= tan(cos(3(1/2)+sin(2*x)./(3+4*x.2); x0=-pi:h:pi;X=-pi:h/12:pi; y0= tan(cos(3(1/2)+sin(2*x0)./(3+4*x0.2); Y= 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);

120、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精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 23 页，共 29 页学习必备欢迎下载); mRXj

121、=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 mRH Rj=x RLj RXj RSj RHj,mRj=mRLj mRXj mRSj mRHj, plot(x0,y0,bo,X,Y,k-,X,YS,rp, X,YH,g), legend(节点 ,被插值函数 ,三次样条函数 ,分段埃尔米特插值函数) 例 6.7.9设函数211)(xxf定义在区间5, 5上，节点 ( X(i),f(X (i) 的横坐标向量 X 的

122、元素是首项a=-5，末项b=5，公差 h=1.5 的等差数列，构造分段三次样条函数)(xSn.把区间5,5分成 20 等份，构成 20 个小区间， 用限定端点约束条件和不限定端点约束条件的MA TLAB 程序计算各小区间中点ix处)(xSn的值，并作出节点， 插值点，)(xf和)(xSn的图形，并与分段埃尔米特插值函数的图形比较. 解 记节点的横坐标,20,2 , 1 , 0, 5. 0,5ihihxi插值点)(21121iiixxx，19,2, 1 ,0i. (1) 不限定端点约束条件，在MATLAB 工作窗口输入程序x0=-5:1.5:5; y0=1./(1+x0.2); x1=-4.75

123、:0.5:4.75; x=-5:0.001:5; y=1./(1+x.2); S= splinetx1(x0,y0,x1,x,y) title( y=1/(1+x2)及其三次样条插值函数, 节点和插值点的图形) 运行后屏幕显示各小区间中点ix处)(xSn的值，作出节点，插值点，)(xf和)(xSn的图形(略). （2）限定端点约束条件，取0)5(1nnSdydy，在 MATLAB工作窗口输入程序x0=-5:1.5:5;y0=1./(1+x0.2);x1=-4.75:0.5:4.75; x=-5:0.001:5; y=1./(1+x.2); S= splinetx (x0,y0,x1,x,y,0

124、,0) title( y=1/(1+x2)及其分段压紧三次样条函数, 节点和插值点的图形) 运行后屏幕显示各小区间中点ix处)(xSn的值（略） ，作出节点， 插值点，)(xf和)(xSn的图形 (略). 如果调节端点约束条件或者增加节点的倍数，例如， 在 MATLAB工作窗口输入程序 x0=-5:0.5:5; y0=1./(1+x0.2); x1=-4.75:0.5:4.75; x=-5:0.001:5; y=1./(1+x.2); S= splinetx (x0,y0,x1,x,y,0,0) title( y=1/(1+x2)及其增加节点后的不限定端点约束条件的三次样条函数, 节点和插值点

125、的图形) 则运行后输出的图形中的三次样条函数与被插值函数的图像基本重合. 例6.7.10设函数xxxfcos5.0)(定义在区间2,2上，取7n，按等距节点构造分段三次样条函数)(xSn， 用MA TLAB 程序计算各小区间中点ix处)(xSn的值，分别作出局部和整体区间上的节点，插值点，)(xf，三次样条函数)(xSn和分段三次埃尔米特插值函数)(xHn的图形，并进行比较. 解编写并保存名为sanci.m 的 M文件如下subplot(2,2,1) h=4*pi/7; x1=-2*pi:h:-2*pi+3*h;c1=0.5.*x1-cos(x1); 精选学习资料 - - - - - - -

126、- - 名师归纳总结 - - - - - - -第 24 页，共 29 页学习必备欢迎下载 xx1 =-2*pi+4*pi/14:h:-2*pi+pi/11+2*h; X1=-2*pi:0.001:-2*pi+3*h; Y1=interp1(x1,c1,xx1, pchip); YY1=interp1(x1,c1,X1, pchip); pp1 = interp1 (x1,c1,xx1,spline ); P1 = interp1 (x1,c1,X1,spline ); cc1 =0.5.*X1-cos(X1);plot(x1,c1,bo,xx1,pp1,k*,X1,P1,k-,xx1,Y1,

127、rx,X1, YY1, r:,X1,cc1,g-.) subplot(2,2,2) x2=2*pi-3*h:h:2*pi;c2=0.5.*x2-cos(x2); xx2 =2*pi-4*pi/14-2*h:h:2*pi-4*pi/14; X2=2*pi-3*h:0.001:2*pi; pp2 = interp1 (x2,c2,xx2, spline ); YY2=interp1(x2,c2, X2, pchip); Y2=interp1(x2,c2,xx2, pchip); P2= interp1 (x2,c2,X2,spline ); cc2=0.5.*X2-cos(X2); plot(x2

128、,c2,bo,xx2,pp2,k*,X2,P2,k-,xx2,Y2,rx,X2, YY2, r:,X2,cc2,g-.) title(y=0.5x-cos(x)及其三次样条函数, 分段三次埃尔米特插值函数,节点和插值点的图形 ) subplot(2,1,2) x=-2*pi:h:2*pi;c=0.5.*x-cos(x); xx =-2*pi+4*pi/14:h:2*pi-4*pi/14; pp = spline(x,c,xx),Y=interp1(x,c,xx, pchip), X=-2*pi:0.001:2*pi; P = interp1 (x,c,X,spline ); YY=interp

129、1(x,c,X, pchip); cc=0.5.*X-cos(X); plot(x,c,bo,xx,pp,k*,X,P,k-,xx,Y,rx,X,YY,r:,X,cc, g-.) legend( 节点 (xi,yi), 三次样条插值 , 三次样条插值函数 , 分段三次埃尔米特插值 , 分段三次埃尔米特插值函数 , y=0.5x-cosx 的函数 ) 在 MATLAB 工作窗口输入文件名 sanci 运行后屏幕显示各小区间中点ix处三次样条插值pp和分段三次埃尔米特插值Y，出现被插值函数)(xf，节点，三次样条和分段三次埃尔米特的函数及其插值点等图形（略）。pp =-3.2181 -0.9609

130、 -0.6824 -0.9476 1.1128 2.6295 2.1675 Y=-3.0085 -1.0074 -0.7589 -0.7872 1.1498 2.4072 2.3787 6.7.8 用 MATLAB 计算有关分段三次样条的误差例6.7.12设函数)432sin3tan(cos()(2xxxf定义在区间,上， 取7n，按等距节点分别作分段线性插值、拉格朗日插值、三次样条插值和分段埃尔米特插值.用MATLAB 程序计算各小区间中点ix处四种插值及其绝对误差和相对误差，作出)(xf及其后三种插值函数，插值点和节点的图形，并进行比较. 解编写并保存名为sanli679.m的 M 文件如

131、下h=2*pi/n; x=(-2*pi+h)/2:h:(2*pi-h)/2; y= tan(cos(sqrt(3)+sin(2*x)./(3+4*x.2); x0=-pi:h:pi;X=-pi:h/12:pi; y0= tan(cos(3(1/2)+sin(2*x0)./(3+4*x0.2); Y= 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); 精选学习资料 - - - - - -

132、- - - 名师归纳总结 - - - - - - -第 25 页，共 29 页学习必备欢迎下载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); mR

133、Xj=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 mRH Rj=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( 节点 , 三次样条插值点 , 被插值函数 , 拉格朗日插值函数 , 三次样条函数 , 分段三次埃尔米特插值函数 ) title(y=

134、tan(cos(sqrt(3)+sin(2x)/(3+4x2)及其三种插值函数 , 节点和插值点的图形 ) 运行程序 n=7;sanli679 得到各小区间中点21ix处的函数值)(21ixy，分段线性插值)(21inxX、拉格朗日插值)(21inxL、三次样条插值)(21inxS和分段埃尔米特插值)(21inxH及其绝对误差、相对误差和平均值的结果, 作出)(xf及其后三种插值函数，插值点和节点的图形(略 ). 例 6.7.13 （机床加工）待加工零件的外形根据工艺要求由一组数据（x, y） 给出（在平面情况下） ，用程控铣床加工时每一刀只能沿x 方向和 y 方向走非常小的一步，这就需要从已

135、知数据得到加工所要求的步长很小的（x, y）坐标 .表 6 15 给出的 x,y 数据位于机翼断面的下轮廓线上（如图 6 25） ，假设需要得到x 坐标每改变0.1 时的 y 坐标 .试完成加工所需数据，画出曲线，并求出x=0 处的曲线斜率和13 x 15 范围内 y 的最小值 . 表 6 15 机翼断面下轮廓线上的部分数据X 0 3 5 7 9 11 12 13 14 15 Y 0 1.2 1.7 2.0 2.1 2.0 1.8 1.2 1.0 1.6 图 6 25机翼断面轮廓线（表6 15 数据用圆点表示）解 根据上述提出的加工要求，以所给数据为节点，在x=0 到 x=15 范围内求步长为

136、0.1 的插值 .用四种插值方法试验，编写并保存名为sancili6710.m程序为 M 文件x0=0 3 5 7 9 11 12 13 14 15 ; x=0:0.1:15; y0=0 1.2 1.7 2.0 2.1 2.0 1.8 1.2 1.0 1.6 ; yL=lagr1(x0,y0,x); xy精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 26 页，共 29 页学习必备欢迎下载yX=interp1(x0,y0,x); yS=interp1(x0,y0,x,spline); yH=interp1(x0,y0,x, pchip); CZ=x

137、 yL yX yS yH subplot(4,1,1) plot(x0,y0,bo,x,yL,r), grid,title( 拉格朗日插值 ) subplot(4,1,2) plot(x0,y0,bo,x,yX,r), grid,title( 分段线性插值 ) subplot(4,1,3) plot(x0,y0,bo,x,yS,r), grid,title( 三次样条 ) subplot(4,1,4) plot(x0,y0,bo,x,yH,r), grid,title( 分段埃尔米特插值 ) 在 MATLAB 工作窗口输入文件名 sancili6710 运行后得到的拉格朗日插值、分段线性插值、

138、三次样条插值和分段埃尔米特插值及其节点的图形，同时还得到拉格朗日插值、分段埃尔米特插值、分段线性插值和三次样条插值的结果 (略). 6.8 高元插值及其 MATLAB 程序6.8.1 MESHGRID 命令的功能和调用格式调用格式一X,Y =meshgrid (x,y)例 6.8.1已知 x=-3:0.2:3；y=x ，计算函数437xze22yx的值，并作出函数的图形 . 解 输入程序 X,Y = meshgrid(-3:.2:3, -3:.2:3); Z =7-3* X.4 .* exp(-X.2 - Y.2), mesh(Z) title( Z =7-3 X4 exp(-X2-Y2)的图

139、形 ) 运行后输出函数值和图形（略）. 例6.8.2作出函数xz2e22yx在区域22x，22y上的图形 .解 输入程序 X,Y = meshgrid(-2:.2:2, -2:.2:2);Z = 2+X .* exp(-X.2 - Y.2); mesh(Z) title(Z = 2+X exp(-X2 - Y2)的图 形) 运行后输出函数值和图形（略）. 调用格式二X,Y = meshgrid (x)调用格式三X,Y,Z = meshgrid (x,y,z) 例6.8.3 计算函数xU2e222zyx在yzyx,5 ,3, 1, 0,4,2,1 ,0处的函数值 . 解输入程序 x=0,1,2,

140、-4;y=0,-1,3,5;z=y; X,Y,Z = meshgrid(x,y,z), U= 2+X .* exp(-X.2 - Y.2- Z.2) 运行后屏幕显示(略).6.8.2 单调数据点上的二元插值及其MATLAB 程序例6.8.4 设节点),(zyx中的值3.0000,1.5000,0,-1.50003.0000-xxy，函数2)1 (3xze)5(1053)1(22yxxyxe3122yxe22)1(yx，作在节点),(zyx处X2,3,1,7,Y5,2,-1,5的双线性插值及其图形. 精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 2

141、7 页，共 29 页学习必备欢迎下载解输入程序 x,y = meshgrid(-3:1.5:3), z=3*(1-x).2.*exp(-(x.2)-(y+1).2)- 10*(x/5 - x.3 - y.5).*exp(-x.2-y.2)-1/3*exp(-(x+1).2-y.2) xi,yi=meshgrid(2,3,1,7,5,2,-1,5); zi=interp2(x,y,z,xi,yi), mesh(xi,yi,zi),xlabel(x), ylabel(y), zlabel(z) title(z=3(1-x)2exp(-x2-(y+1)2)-10(x/5-x3-y5)exp(-x2-

142、y2)- 1/3 exp(-(x+1)2-y2) 的双线性插值图形) 运行后屏幕显示双线性插值及其图形( 略) . 例6.8.5设节点),(zyx中的3:0.5:3-x,xy和函数337xZe22yx，作Z在插值点 X,5:5.0:.93-Y5.4:5. 0:9 .4-处的二元样条插值， 双三次插值和数据点的图形 .解（1）计算二元样条插值. 输入程序 x,y = meshgrid(-3:0.5:3); z =7-3* x.3 .* exp(-x.2 - y.2); xi=-3.9:0.5:5;yi=-4.9:0.5:4.5; xi,yi = meshgrid(xi,yi); zi = int

143、erp2(x,y,z,xi,yi, spline), mesh(xi,yi,zi), hold on, plot3(x,y,z,r. ,markersize,3*5), hold off xlabel(x ), ylabel(y ), title(z =7-3 x3 exp(-x2 - y2) 的二元样条插和值数据点的图形) 运行后屏幕显示Z在插值点 X,5:5.0:.93-Y 5. 4:5.0:9. 4-处的二元样条插值及其图形 ( 略) .（2）计算双三次插值. 输入程序 x,y = meshgrid(-3:0.5:3); z =7-3* x.3 .* exp(-x.2 - y.2); x

144、i=-3.9:0.5:5; yi=-4.9:0.5:4.5; xi,yi = meshgrid(xi,yi); zi=interp2(x,y,z,xi,yi, cubic), mesh(xi,yi,zi),hold on plot3(x,y,z,r. ,markersize,3*5), hold off xlabel(x), ylabel(y), zlabel(z), title(z=7-3x3exp(-x2-y2) 的双三次插值和数据点(x,y,z)的图形 ) 运行后屏幕显示Z在插值点 X,5:5.0:.93-Y5. 4:5.0:9. 4-处的双三次插值和数据点图形（略） （三种方法比较留给

145、读者）例 6.8.6 设节点),(zyx中的5:0.5:5-x,xy和函数437xZe22yx，作Z在插值点X,5:5.0:.93-Y5.4:5.0:9 .4-处的双三次插值和二元最近邻插值及其图形 . 解（1）双三次插值.输入程序 x,y = meshgrid(-5:0.5:5); z =7-3* x.4 .* exp(-x.2 - y.2); xi=-3.9:0.5:5; yi=-4.9:0.5:4.5; xi,yi = meshgrid(xi,yi); zi = interp2(x,y,z,xi,yi, cubic), mesh(xi,yi,zi) hold on plot3(x,y,z

146、,r., markersize,3*5) hold off xlabel(x), ylabel(y), zlabel(z), title(z =7-3 x4 exp(-x2 - y2) 的双三次插值和数据点的图形 ) 运行后屏幕显示Z在插值点X, 5:5 .0:.93-Y 5. 4:5 .0:9. 4-处的双三次插值及其图形(略). 精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 28 页，共 29 页学习必备欢迎下载（2）二元最近邻插值. 输入程序 x,y = meshgrid(-5:0.5:5); z =7-3* x.4 .* exp(-x.2

147、 - y.2); xi=-3.9:0.5:5; yi=-4.9:0.5:4.5; xi,yi = meshgrid(xi,yi); zi = interp2(x,y,z,xi,yi, nearest), mesh(xi,yi,zi) hold on,plot3(x,y,z,r., markersize,3*5), hold off xlabel(x), ylabel(y), zlabel(z) title(z =7-3 x3 exp(-x2-y2) 的二元最近邻插值和数据点的图形 ) 运行后屏幕显示Z在插值点X, 5:5.0:.93-Y5. 4:5.0:9. 4-处的二元最近邻插值及其图形 (

148、 略). 6.8.3 三元插值及其 MATLAB 程序例6.8.7 设节点),(zyx的坐标为yzyx,15,3 ,0 ,1,12, 1 , 0, 4，计算函数xV2e222zyx在插值点13:25.0:3,3:25.0:3,10:25.0:3iiizyx处的三元线性插值，并作其图形. 解输入程序 x=-4,0,1,12;y=-1,0,3,15;z=y; X,Y,Z= meshgrid(x,y,z); V= 2+X .* exp(-X.2 - Y.2- Z.2); xi,yi,zi = meshgrid(-3:.25:10,-3:.25:3,-3:.25:13); vi = interp3(X

149、,Y,Z,V,xi,yi,zi), slice(xi,yi,zi,vi,-1 6 9.5,9,-2 .2 9), shading flat,lighting flat xlabel(x), ylabel(y), zlabel(z), title(V=2+ x exp(-x2 - y2-z2) 的三元线性插值图形 ) hold on,colorbar(horiz), view(-30 45) 运行后屏幕显示三元线性插值及其图形（略）.例6.8.8 取n=10 ，作函数 flow 在插值点,3:25.0:3,10:25.0:1 .0iiyxiiyz处的三元三次样条插值及其图形. 解输入程序 x,y

150、,z,v = flow(10); xi,yi,zi = meshgrid(.1:.25:10,-3:.25:3,-3:.25:3); vi = interp3(x,y,z,v,xi,yi,zi,spline); slice(xi,yi,zi,vi,2.5 7.5,0.1,-1.2 1.5), shading flat,xlabel(x), ylabel(y), zlabel(z), title( flow 的三元三次样条插值图形 ) hold on,colorbar(horiz) 运行后屏幕显示三元三次样条插值及其图形（略）.精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 29 页，共 29 页