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

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

1、学习必备欢迎下载2.1 插值问题及其误差2.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)例 2.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) 2.2 拉格朗日（ Lagrange）插值及其 MATLAB 程序2.2.1 线性插值及其 MATLAB 程序例 2.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为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 运行后输出误差限为 R1 = 精选学习资料 - - - -

5、 - - - - - 名师归纳总结 - - - - - - -第 1 页，共 21 页学习必备欢迎下载 1.8750 例 2.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 和插值多项式L 为l

6、0 = 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. 2.2.2 抛物线插值及其MATLAB 程序例2.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 = 5734161139222659/90

8、07199254740992*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) , 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

9、(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), 精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 2 页，共 21 页学习必备欢迎下载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= -6048313895780875/18014

11、398509481984*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. 2.2.3 n次拉格朗日（ Lagrange）插值及其 MATLAB 程序例2.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,1,2,17;

12、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 ( conv (p1,

13、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*x3+1/8*x2-

14、1/4*x 输入程序 L=poly2sym (P),x=0.6; Y = polyval(P,x) 运行后输出插值多项式和插值为L = 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 运行后输出误差限为R3 = 91/2500*M 即精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 3 页，共 21 页学习必备欢迎下载R3)(04.0)4(f，)00.2,00.2(. 2.2.5 拉格朗日多项式和基函数的MATLAB

15、程序求拉格朗日插值多项式和基函数的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例2.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(f， 作五次拉格

16、朗日插值多项式和基函数，并写出估计其误差的公式.解在 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 = -0.005

17、6 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.0323*x3-0

18、.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+0.0220*x

19、-0.0002 估计其误差的公式为)(5xR)25. 3)(03. 2)(02. 1)(01.0()00. 1)(15. 2(! 6)()6(xxxxxxf，)3.25,-2.15(.2.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; for k=1:n 精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 4 页，共 21 页学习必备欢迎下载

20、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; 例 2.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.8660; y,R

21、=lagranzi(X,Y,x,M) 运行后输出插值y 及其误差限R 为y = R = 0.6434 8.8610e-004. 2.53 分段插值及其 MATLAB 程序2.3.1 高次插值的振荡和MATLAB 程序例 2.3.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

22、.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: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

23、); plot(x,y,k-), gtext(y=tan(cos(sqrt(3)+sin(2x)/(3+4x2),pause n=3; x0=-pi:2*pi/(3-1):pi; 精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 5 页，共 21 页学习必备欢迎下载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的图形 (略).

24、 （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./(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=

25、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: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

26、(x,y5,m ),gtext(n=10) title( 高次拉格朗日插值的振荡 ) 回车运行后，便会逐次画出)1/(12xy在区间5,5上的n次拉格朗日插值多项式)(xLn)10,8,6,4,2(n的图形 (略).2.3.4 作有关分段线性插值图形的MATLAB 程序作有关分段线性插值图形的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), 分段线性插值函数

27、Sn(x), 插值点 (x,s), 精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 6 页，共 21 页学习必备欢迎下载被插值函数 y ) 我们也可以直接在在MATLAB 工作窗口编程序. 例如， x0 =-6:6; y0 =sin(x0); xi = -6:.25:6; yi = interp1(x0,y0,xi); x=-6:0.001:6; y=sin(x);plot(x0,y0,o,xi,yi,x,y,:), legend( 节点 (xi,yi), 分段线性插值函数 , 被插值函数 y=sinx ) title(y=sinx及其分段线性插值

28、函数和节点的图形) x0 =-6:6; y0 =cos(x0); xi = -6:.25:6; yi = interp1(x0,y0,xi); x=-6:0.001:6; y=cos(x); plot(x0,y0,o,xi,yi,x,y,:), legend( 节点 (xi,yi), 分段线性插值函数 , 被插值函数 y=cosx ) title(y=cosx及其分段线性插值函数和节点的图形)例2.3.5设函数22511)(xxf， 在区间 1 , 1上取等距节点),(iiyx10, 2, 1, 0i, 构造分段线性插值函数)(xSn，用 MATLAB 程序计算各小区间的中点ix处)(xSn的

29、值，作出节点，插值点，)(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=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 throug

30、h 4 0.04864253393665 0.07941176470588 0.15000000000000 0.35000000000000 Columns 5 through 8 0.75000000000000 0.75000000000000 0.35000000000000 0.15000000000000 Columns 9 through 10 0.07941176470588 0.04864253393665 2.3.5 用 MATLAB 计算有关分段线性插值的误差例2.3.8设函数)432sin3tan(cos()(2xxxf,在区间,上取等距节点),(iiyx，9, 2,

31、1, 0i，构造分段线性插值函数)(xSn. (1)用MATLAB 程序计算各小区间中点ix处)(xSn的值，作出节点， 插值点，)(xf和)(xSn的图形；(2) 用MATLAB 程序计算各小区间中点处)(xSn的值及其相对误差；(3) 用MATLAB 程序估计)(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; 精选学习资料 - - - - - - - - - 名师归纳总结 - - -

32、- - - -第 7 页，共 21 页学习必备欢迎下载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,Ri,i=xi,fi,si,Ri title(y=tan(cos(sqrt(3)+sin(2 x)/(3+4 x2)的分段线性插值的有关图形 ) 运行后屏幕显示Ri（略），并且作出节点，插值点，)(xf和

33、)(xSn的图形 (略). 在 MATLAB 工作窗口输入程序 syms x y=tan(cos(3(1/2)+sin(2*x)/(3+4*x2);yxx=diff(y,x,2), 运行后屏幕显示函数)(xf的二阶导数)( xf（略） . 在 MATLAB 工作窗口输入程序 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.

34、*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).*(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.*

35、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)* myxx/8 运行后屏幕显示myxx =)(maxxfx和)(xSn在区间,上的误差限R如下myxx = R = -0.02788637150664 -0.00169893490711 2.4 分段埃尔米特插值及其MATLAB 程序2.4.2 分段埃尔米特插值的MATLAB 程序调用格式一：YI=interp1(X,Y,XI,p

36、chip) 调用格式二 ：YI=interp1(X,XI,pchip)例 2.4.5试用 MA TLAB 程序计算例6.6.1中在各小区间中点处分段三次埃尔米特插值)(2/1inxH及其相对误差 . 解在 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); yi=interp1(x0,y0,xi,pchip); Ri=abs(fi-yi)./fi); xi,fi,yi,Ri,i=xi,fi,yi,Ri 运行后屏幕显示各小区间中点xi处的函数值fi，插值 si，相对误差值Ri

37、 (略). 2.4.3 作有关分段埃尔米特插值图形的MATLAB 程序作有关分段埃尔米特插值图形的MATLAB 主程序function H=hermitetx(x0,y0,xi,x,y) H= interp1(x0,y0,xi,pchip); Hn= interp1(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

38、:6; 精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 8 页，共 21 页学习必备欢迎下载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) title( y=sinx及其分段埃尔米特插值函数和节点的图形 ) x0 =-6:6; y0 =cos(x0); xi = -6:.25:6;yi = interp1(x0,y0,xi,pchip); x=-

39、6:0.001:6; y=cos(x); plot(x0,y0,o,xi,yi,x,y,:), legend( 节点 (xi,yi), 分段埃尔米特插值函数 , 被插值函数y=cosx) title( y=cosx及其分段埃尔米特插值函数和节点的图形) 例2.4.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的值，并作出节点，插值点，

40、)(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,x1,x,y) title( 函数 y=1/(1+x2)及其分段埃尔米特插值函数，插值，节点 (xi,yi) 的图形 )运行后屏幕显示各小区间中点ix处)(3,xHn的值，出现节点，插值点，)(xf和)(3,xHn的图形 (略).例2.4.7设函数xxxfcos5.0)(定义在区间,上，取7n，按等距节点构造分段埃尔米特插值函数)(3,7

41、xH，用 MATLAB 程序计算各小区间中点ix处)(3,7xH的值，作出节点，插值点，)(xf和)(3, 7xH的图形 .解记节点的横 坐标7,2 ,1 ,0,7/2,ihihxi插 值点)(21121iiixxx，6,2, 1 ,0i.在 MATLAB工作窗口输入程序 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(

42、x)及其分段埃尔米特插值函数，插值，节点(xi,yi) 的图形 )运行后屏幕显示各小区间中点ix处)(3,7xH的值，出现节点， 插值点，)(xf和)(3,7xH的图形 (略).2.4.4 用 MATLAB 计算有关分段埃尔米特插值的误差例2.4.8 设函数22511)(xxf定义在区间 1,1上，取10n，按等距节点构造分段埃尔米特插值函数)(3,xHn，用 MATLAB 程序在 1,1上计算)(max)4(11xfx和)(3,xHn的误差公式和误差限.解在 MATLAB 工作窗口输入程序精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 9 页，共

43、 21 页学习必备欢迎下载 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上绝对值的最大值myxxxx 和)(3,xHn在区间 1, 1上的误差公式myxxxx 为myxxxx = R = 15000 625/16*h4 （4）在 MATLAB 工作窗口输入程序 h=0.2; R =625/16

44、*h4 运行后输出误差限为R = 0.06250000000000 例2.4.9 设函数)432sin3tan(cos()(2xxxf定义在区间,上，取9n，按等距节点构造分段埃尔米特插值函数)(3,xHn. (1)用MATLAB 程序计算各小区间中点ix处)(3,xHn的值， 作出节点， 插值点，)(xf和)(3,xHn的图形；(2) 并用 MATLAB 程序计算各小区间中点处)(3,xHn的值及其相对误差；(3) 用MATLAB 程序求)(max)4(xfx和)(3,xHn在区间,上的误差公式和各插值的误差限.解 （1）记节点的横坐标9 , 2, 1 , 0, 9/2,ihihxi，插值点

45、)(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; 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,f

46、i,Hi,Ri,i=xi,fi,Hi,Ri title( 函数 y= tan(cos(sqrt(3)+sin(2x)/(3+4x2)及其分段埃尔米特插值函数，插值，节点(xi,yi) 的图形 )运行后屏幕显示各小区间中点xi处的函数值fi，插值 Hi，相对误差值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，

47、 然 后 将 输 出 的)()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.*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

48、.*x)./(3.+4.*x.2).2.*x+128.*(3.(1./2)+sin(2.*x)./(3.+4.*x.2).3.*x.2.-8.精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 10 页，共 21 页学习必备欢迎下载*(3.(1./2)+sin(2.*x)./(3.+4.*x.2).2)+16.*(1.+tan(cos(3.(1./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.

49、(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)./(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)

50、+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).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.

51、/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)+(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(

52、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+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)

53、.*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.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+

54、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.*cos(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

55、./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(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

56、)./(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)./(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

57、+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)./(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

58、(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)./(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.

59、*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.*x)./(3+4.*x.2).2).*cos(3.(1./2)+sin(2.*x)./(3精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 11 页，共 21 页学习必备欢迎下载+4.*x.2).2.*(2.*cos(2.*x)./(3+4.*x.2)-8.*(3.(1./2)+sin(2.*

60、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.*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(c

61、os(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(2.*x)./(3+4.*x.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.*

62、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./384) 将其保存为名为myxx.m的 M文件，然后在MATLAB 工作窗口输入该文件名 myxx运 行 后 屏 幕 显 示myxx =)(max)4(xfx和)(3,xHn在 区 间,上 的 误 差 公 式)(ma x384)4(4xfhRx如下myxx = R = 73.94706841647552 173452078002906

63、1/9007199254740992*h4 最后在 MATLAB 工作窗口输入h=2*pi/9; R =1734520780029061/9007199254740992*h4 运行后屏幕显示)(3,xHn在区间,上的误差限R = 0.04574453029948 2.5 三次样条（ Spline ）及其 MATLAB 程序2.5.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

64、=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);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

65、+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,m)=2; A(i,i+1)=lambda(i); A(i+1,i)=mu(i); end dY=AD; 精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 12 页，共 21 页学习必备欢迎下载syms x m=length(X);S=zeros(m-1,1); for k=2:m sk=Y(k

66、-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例 2.5.2 用计算压紧样条的MATLAB 主程序 ClampedSp.m求例 6.7.1的问题 . 解保存 ClampedSp.m为 M 文件，输入程序 X=0:2:6,10,Y=0,16,36,54,82,dyo=8,dyn=7, m,H,lambda,

67、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，线性方程组（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*

68、(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-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

69、 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)*(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

70、()8(4sf的值如下精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 13 页，共 21 页学习必备欢迎下载x2 = s2 = x4 = s4 = 3 25.7500 8 68.5000 例 2.5.2和例 2.5.1计算的结果完全相同. 计算 自然样条的MATLAB 主程序function m,H,lambda,mu,D,A,dY,sk=NaturalSp(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

71、(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; 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)=d

72、k; 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); 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)

73、*(x-X(k)2)/(H(k-1)2)+dY(k)*(x-X(k)*(x-X(k-1)2)/(H(k-1)2) end 例 2.5.3 用计算自然样条的MATLAB 主程序 NaturalSp.m求例 6.7.1的问题 . 解保存 ClampedSp.m为 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， （

74、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/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)*(

75、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 A = 精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 14 页，共 21 页学习必备欢迎下载 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

76、/43 输入程序 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+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 例 2.5.3和例 2.5.2用的方法不同，所以计算的结果不同，但是

77、两种方法计算的结果相近 . 2.5.6 用 MATLAB 计算三次样条例 2.5.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.65 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

78、.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), 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.

79、5000 coefs: 8x4 double pieces: 8 精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 15 页，共 21 页学习必备欢迎下载order: 4 dim: 1（2）因为端点约束条件为0) 1(nS，0)50.3(nS，所以输入MATLAB 程序 YI= spline (X, 0,Y,0,XI), PP= spline (X, 0,Y,0) 运行后屏幕显示压紧三次样条分别在02.01x，56. 22x的插值和该样条函数的有关信息如下YI = -3.0192 4.7501 PP = form: pp breaks: -1 -0.

80、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.5.7在默认端点约束条件下，用两种方法计算在下列插值点处的三次样条插值. （1）给出节点的数据与例2.5.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个等分点 .

81、解 （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= 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 = l

82、inspace(-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.5935 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 例 2.5.8设函数22511)(xxf定义在区间 1, 1上，取10n，按等距节点构造分段三次样条函数)(xSn.试用MAT

83、LAB程序计算在各小区间中点处分段三次样条插值)(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); yi=spline (x0,y0,xi); Ri=abs(fi-yi)./fi); xi,fi,yi,Ri,i=xi,fi,yi,Ri 运行后屏幕显示各小区间中点xi处的函数值fi，插值 si，相对误差值Ri(略). 2.5.7 几种作三次样条有关图像的MATLAB 程序精选学习资料 - - - - - - - - - 名师归纳总结 -

84、- - - - - -第 16 页，共 21 页学习必备欢迎下载求有关分段三次样条图形的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, *,x,y,-.) legend( 节点 (xi,yi), 分段三次样条函数 , 插值点 (x,S), 被插值函数 y ) （二）不限定端点约束条件的作图程序function S=splinetx1(x0,

85、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), 分段三次样条函数 , 插值点 (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 =

86、 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 (x2,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 = spli

87、ne(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,c1=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)

88、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,bo,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

89、,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,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); 精选学习资料 - - - - - - - -

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

91、./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 mRH Rj=x RLj RXj RSj RHj,mRj=mRLj mRXj mRSj mRHj, plo

92、t(x0,y0,bo,X,Y,k-,X,YS,rp, X,YH,g), legend(节点 ,被插值函数 ,三次样条函数 ,分段埃尔米特插值函数) 例 2.5.9设函数211)(xxf定义在区间5, 5上，节点 ( X(i),f(X (i) 的横坐标向量 X 的元素是首项a=-5，末项b=5，公差 h=1.5 的等差数列，构造分段三次样条函数)(xSn.把区间5,5分成 20 等份，构成 20 个小区间， 用限定端点约束条件和不限定端点约束条件的MA TLAB 程序计算各小区间中点ix处)(xSn的值，并作出节点， 插值点，)(xf和)(xSn的图形，并与分段埃尔米特插值函数的图形比较. 解

93、记节点的横坐标,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: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）限定端点约束条件，

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

95、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)及其增加节点后的不限定端点约束条件的三次样条函数, 节点和插值点的图形) 则运行后输出的图形中的三次样条函数与被插值函数的图像基本重合. 精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 18 页，共 21 页学习必备欢迎下载例2.5.10设函数xxxfcos5.0)(定义在区间2,2上，取7n，按等距节点构造分段三次样条函数)(xSn， 用MA TLAB

96、 程序计算各小区间中点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); 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);

97、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,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, s

98、pline ); 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,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*p

99、i+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=interp1(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

100、的函数 ) 在 MATLAB 工作窗口输入文件名 sanci 运行后屏幕显示各小区间中点ix处三次样条插值pp和分段三次埃尔米特插值Y，出现被插值函数)(xf，节点，三次样条和分段三次埃尔米特的函数及其插值点等图形（略）。pp =-3.2181 -0.9609 -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 2.5.8 用 MATLAB 计算有关分段三次样条的误差例2.5.12设函数)432sin3tan(cos()(2xxxf定义在区间,上， 取7n，按等距

101、节点分别作分段线性插值、拉格朗日插值、三次样条插值和分段埃尔米特插值.用精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 19 页，共 21 页学习必备欢迎下载MATLAB 程序计算各小区间中点ix处四种插值及其绝对误差和相对误差，作出)(xf及其后三种插值函数，插值点和节点的图形，并进行比较. 解编写并保存名为sanli679.m的 M 文件如下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

102、; 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); 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

103、)./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 mRH Rj=x RLj RXj RSj RHj,mRj=mRLj mRXj mRSj mRHj, pl

104、ot(x0,y0,bo,x,yS,r*,X,Y,k-,X,YL,g-.,X,YS,c:, X,YH, m-), legend( 节点 , 三次样条插值点 , 被插值函数 , 拉格朗日插值函数 , 三次样条函数 , 分段三次埃尔米特插值函数 ) title(y=tan(cos(sqrt(3)+sin(2x)/(3+4x2)及其三种插值函数 , 节点和插值点的图形 ) 运行程序 n=7;sanli679 得到各小区间中点21ix处的函数值)(21ixy，分段线性插值)(21inxX、拉格朗日插值)(21inxL、三次样条插值)(21inxS和分段埃尔米特插值)(21inxH及其绝对误差、相对误差和

105、平均值的结果, 作出)(xf及其后三种插值函数，插值点和节点的图形(略 ). 例 2.5.13 （机床加工）待加工零件的外形根据工艺要求由一组数据（x, y） 给出（在平面情况下） ，用程控铣床加工时每一刀只能沿x 方向和 y 方向走非常小的一步，这就需要从已知数据得到加工所要求的步长很小的（x, y）坐标 .表 6 15 给出的 x,y 数据位于机翼断面的下轮廓线上（如图 6 25） ，假设需要得到x 坐标每改变0.1 时的 y 坐标 .试完成加工所需数据，画出曲线，并求出x=0 处的曲线斜率和13 x 15 范围内 y 的最小值 . 机翼断面下轮廓线上的部分数据X 0 3 5 7 9 11

106、 12 13 14 15 Y 0 1.2 1.7 2.0 2.1 2.0 1.8 1.2 1.0 1.6 xy精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 20 页，共 21 页学习必备欢迎下载机翼断面轮廓线解 根据上述提出的加工要求，以所给数据为节点，在x=0 到 x=15 范围内求步长为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

107、; yL=lagr1(x0,y0,x); yX=interp1(x0,y0,x); yS=interp1(x0,y0,x,spline); yH=interp1(x0,y0,x, pchip); CZ=x 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 运行后得到的拉格朗日插值、分段线性插值、三次样条插值和分段埃尔米特插值及其节点的图形，同时还得到拉格朗日插值、分段埃尔米特插值、分段线性插值和三次样条插值的结果 (略). 精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 21 页，共 21 页