第六章 函数的插值方法 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 求三个一次多项式g(x)、h(x)和f(x)的积f(x)g(x)h(x).它们的零点分别依次为0.4,0.8,1.2.
解 我们可以用两种MATLAB程序求之. 方法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 =
x^3-12/5*x^2+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] =deconv (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 已知函数f(x)在[1,3]上具有二阶连续导数，f\"(x)5，且满足条件f(1)1,f(3)2.求线性插值多项式和函数值f(1.5)，并估计其误差.
解 输入程序
>> X=[1,3];Y=[1,2]; 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=1.5; Y = polyval(P,x)
运行后输出基函数l0和l1及其插值多项式的系数向量P（略）、插值多项式L和插值Y为
76．
第三篇 第六章 函数的插值方法
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 =
1.8750
例6.2.2 求函数f(x)e在[0,1]上线性插值多项式，并估计其误差. 解 输入程序
>> 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为
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.
x
6.2.2 抛物线插值及其MATLAB程序
例6.2.3 求将区间 [0, π/2] 分成n等份(n1,2)，用yf(x)cosx产生
n1个节点，然后根据（6.9）和（6.13）式分别作线性插值函数P1(x)和抛物线插值函数P2(x).用它们分别计算cos (π/6) (取四位有效数字)，并估计其误差.
解 输入程序
>> 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/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) ,
77．
第三篇 第六章 函数的插值方法
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), P = l01* Y(1)+ l11* Y(2) + l21* Y(3), L=poly2sym (P),x=pi/6; Y = polyval(P,x)
运行后输出基函数l01、l11和l21及其插值多项式的系数向量P、插值多项式L和插值Y为
l0 =
228155022448185/281474976710656*x^2-2150310427208497/1125899906842624*x+1
l1 =
-228155022448185/140737488355328*x^2+5734161139222659/2251799813685248*x
l2 =
228155022448185/281474976710656*x^2-5734161139222659/9007199254740992*x
P =
-0.3357 -0.1092 1.0000 L=
-6048313895780875/18014398509481984*x^2-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 给出节点数据f(2.00)17.00，f(0.00)1.00，f(1.00)2.00，
f(2.00)17.00，作三次拉格朗日插值多项式计算f(0.6)，并估计其误差.
解 输入程序
>> 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(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),
运行后输出基函数l0，l1，l2和l3及其插值多项式的系数向量P（略）为
l0 =
-1/24*x^3+1/8*x^2-1/12*x，l1 =1/4*x^3-1/4*x^2-x+1 l2 =
-1/3*x^3+4/3*x，l3 =1/8*x^3+1/8*x^2-1/4*x
输入程序
>> L=poly2sym (P),x=0.6; Y = polyval(P,x)
运行后输出插值多项式和插值为
L = Y =
x^3+4*x^2-4*x+1 0.2560.
78．
第三篇 第六章 函数的插值方法
输入程序
>> 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 0.04f(4)()， (2.00,2.00).
6.2.5 拉格朗日多项式和基函数的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 给出节点数据f(2.15)17.03，f(1.00)7.24，f(0.01)1.05，
f(2.03)17.06，f(3.25)23.05，作五次拉格朗日插值多项式和基函数，f(1.02)2.03，
并写出估计其误差的公式.
解 在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*x^2+2.1076*x^3+0.0648*x^4-0.2169*x^5 L1 =
-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*x^5+0.0299*x^4-0.0323*x^3-0.0292*x^2+0.0382*x-0.0004] [ 0.0331*x^5-0.1377*x^4-0.0503*x^3+0.6305*x^2-0.4852*x+0.0048] [ -0.0693*x^5+0.2184*x^4+0.3961*x^3-1.2116*x^2-0.3166*x+1.0033] [ 0.0687*x^5-0.1469*x^4-0.5398*x^3+0.6528*x^2+0.9673*x-0.0097] [ -0.0317*x^5+0.0358*x^4+0.2530*x^3-0.0426*x^2-0.2257*x+0.0023] [ 0.0049*x^5+0.0004 *x^4-0.0266*x^3+0.0001*x^2+0.0220*x-0.0002]
估计其误差的公式为
f(6)()R5(x),3.25). (x2.15)(x1.00)(x0.01)(x1.02)(x2.03)(x3.25)，(-2.156!
79．
第三篇 第六章 函数的插值方法
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; for 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 已知sin300.5，sin450.7071，sin600.8660，用拉格朗日插值及其误差估计的MATLAB主程序求sin40的近似值，并估计其误差.
解 在MATLAB工作窗口输入程序
>> x=2*pi/9; M=1; X=[pi/6 ,pi/4, pi/3];
Y=[0.5,0.7071,0.8660]; [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,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 给出节点数据f(2.15)17.03，f(1.00)7.24，f(0.01)1.05，f(1.02)2.03，f(2.03)17.06，f(3.25)23.05,作五阶牛顿插值多项式和差商，并写出其估计误差的公式.
解 （1）保存名为newpoly.m的M文件. （2）输入MATLAB程序
80．
第三篇 第六章 函数的插值方法
>> 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及其系数向量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/36028797018963968*x^5+583849564517807/9007199
254740992*x^4+593245028711263/281474976710656*x^3+3823593773002357/1125899906842624*x^2-321902673270315/70368744177664*x+308328649211299/281474976710656
wcgs =
1/720*M*(x^6-79/25*x^5-14201/2500*x^4+4934097026900981/2814749767
10656*x^3+154500712237335/35184372088832*x^2-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*x^2+2.1076*x^3+0.0648*x^4-0.2169*x^5.
估计其误差的公式为
f(6)(),3.25). (x2.15)(x1.00)(x0.01)(x1.02)(x2.03)(x3.25)，(-2.15R5(x)6!
例6.3.4 求函数f(x)7e
x/5在[2,6]上五阶牛顿插值多项式，估计其误差的公
式和误差限公式.用它们计算f(3.156)，并估计其误差.
解 （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)/(5^6)),
运行后输出差商矩阵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.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.0000 C =
0.0000 -0.0004 0.0089 -0.1389 1.3985 -6.9991 L =
9721799720875/1152921504606846976*x^5-3503994098647815/9223
372036854775808*x^4+160742008798419/18014398509481984*x^3-1251152213853501/9007199254740992*x^2+6298131904328647/4503599627370496*x-3940156929554013/562949953421312
wcgs =
1/720*M*(x^6-24*x^5+1172/5*x^4-5952/5*x^3+7276634802928539/219
9023255552*x^2-5237461186650519/1099511627776*x+6085939356447121/219902
81．
第三篇 第六章 函数的插值方法
3255552)
Cw =
0.0014 -0.0333 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*(x^6-24*x^5+1172/5*x^4-5952/5*x^3
+7276634802928539/2199023255552*x^2-5237461186650519/1099511627776*x+6085939356447121/2199023255552), 运行后输出误差限公式wcgs1如下
wcgs1 =
5565367633581443/158456325028528675187087900672*x^6-16696102900744329/19807040628566084398385987584*x^5+1630652716639362799/198070406285660843983859875840*x^4-517579189923074199/12379400392853802748991242240*x^3+40497147813610772932297365501777/348449143727040986586495598010130648530944*x^2-29148396970323855270001164718917/174224571863520493293247799005065324265472*x+33870489914310283926665276375603/348449143727040986586495598010130648530944
（3）输入MATLAB程序
>> x=3.156;
y=9721799720875/1152921504606846976*x^5-3503994098647815
/9223372036854775808*x^4+160742008798419/18014398509481984*x^3-1251152213853501/9007199254740992*x^2+6298131904328647/4503599627370496*x-3940156929554013/562949953421312
wcgs2=1/720*M*(x^6-24*x^5+1172/5*x^4-5952/5*x^3+72766348
02928539/2199023255552*x^2-5237461186650519/1099511627776*x+6085939356447121/2199023255552)
运行后输出f(3.156)的近似值y，及其误差限wcgs2如下
y =
-3.7237 wcgs2 =
-2.4764e-007
6.3.4 牛顿插值及其误差估计的MATLAB程序
牛顿插值及其误差估计的MATLAB主程序
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=length(C); C(d)=C(d)+A(k,k); end
y(k)= polyval(C, z); end
R=M*q1/c1;
例6.3.5 已知sin300.5，sin450.7071，sin600.8660，用牛顿插值
82．
第三篇 第六章 函数的插值方法
法求sin40的近似值，估计其误差，并与例 6.2.6的计算结果比较.
解 方法1（牛顿插值及其误差估计的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 = 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 =
-1583578379808357/4503599627370496*x^2+1408883391907005/
1125899906842624*x-132405829044691/2251799813685248
wcgs =
1/6*M*(x^3-3/4*x^2*pi+4012734077357799/2251799813685248*
x-7757769783530263/18014398509481984)
Cw =
0.1667 -0.3927 0.2970 -0.0718 y =
0.6434
上述两种方法计算y的结果相同.
6.3.5 牛顿插值法的MATLAB综合程序
求牛顿插值多项式、差商、插值及其误差估计的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=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 给出节点数据f(4.00)27.00，f(0.00)1.00，f(1.00)2.00，
f(2.00)17.00，作三阶牛顿插值多项式，计算f(2.345)，并估计其误差.
解 首先将名为newdscg.m的程序保存为M文件，然后在MATLAB工作窗口输入程序
>> 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)
83．
第三篇 第六章 函数的插值方法
运行后输出插值yf(2.345)及其误差限公式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.9167 4.2500 -4.1667 1.0000 P =
11/12*x^3+17/4*x^2-25/6*x+1
例6.3.7 将区间 [0,π/2] 分成n等份(n2,3)，用yf(x)sinx产生n1个节点，求二阶和三阶牛顿插值多项式P2(x)和P3(x).用它们分别计算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]=newdscg(X1,Y1,x,M),
[y2,R2,A2,C2,P3]=newdscg(X2,Y2,x,M)
运行后输出插值y1=P2(/7)sin(/7)和y2=P3(/7)sin(/7)及其误差限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*x^2+163820246322191/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*x^3-4717668559161127/7
2057594037927936*x^2+4595602396951547/4503599627370496*x
6.4 埃尔米特（Hermite）插值及其MATLAB程序
84．
第三篇 第六章 函数的插值方法
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(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 给定函数f(x)在点x0/6,x1/4,x2/2处的函数值f(x0)0.5,
f(x1)0.7071，
f(x2)1和导数值f'(x0)0.8660，f'(x1)0.7071，
(6)f'(x0)0.0000，且f(x)1，求函数f(x)在点x0,x1,x2处的五阶埃尔米特插值
)并估计其误差. 多项式H5(x)和误差公式，计算f(1.567解 （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 =
1.0e+003 *
0.1912 -0.9273 1.6903 -1.4380 0.5751 -0.0866 Hk =
6725828781679091/35184372088832*x^5-4078286086775209/4398
046511104*x^4+7434035571017927/4398046511104*x^3-3162205449085973/2199023255552*x^2+5058863928652835/8796093022208*x-6094057839958843/70368744177664
wcgs =
1/720*M*(x^6-11/6*x^5*pi+7446708432019761/562949953421312
*x^4-4363745503235773/281474976710656*x^3+21569239021155/2199023255552*x^2-7178073637328281/2251799813685248*x+3758430567659515/9007199254740992)
Cw =
0.0014 -0.0080 0.0184 -0.0215 0.0136 -0.0044 0.0006 当f(6)(x)1M时的误差公式为
R=0.001 4* x^6-0.008 0*x^5+0.018 4*x^4-0.021 5*x^3+0.013 6*x^2-0.004 4*x+0.000 6
（3）在MATLAB工作窗口输入程序
>> x=1.567;M=1;
Hk=6725828781679091/35184372088832*x^5-4078286086775209/
4398046511104*x^4+7434035571017927/4398046511104*x^3-3162205449085973/2199023255552*x^2+5058863928652835/8796093022208*x-6094057
85．
第三篇 第六章 函数的插值方法
839958843/70368744177664, wcgs=1/720*M*(x^6-11/6*x^5*pi+7446708432019761/562949953
421312*x^4-4363745503235773/281474976710656*x^3+21569239021155/2199023255552*x^2-7178073637328281/2251799813685248*x+3758430567659515/9007199254740992)
)的近似值Hk及其误差wcgs如下 运行后输出f(1.567Hk =
2.5265 wcgs =
1.3313e-008
6.5 分段插值及其MATLAB程序
6.5.1 高次插值的振荡和MATLAB程序
例6.5.1 作下列函数在指定区间上的n次拉格朗日插值多项式Ln(x) (n2,4,6,8,10)的图形，并讨论插值多项式Ln(x)的次数与误差Rn(x)的关系.
（1）f(x)tan(cos(3sin2x1f(x)x[,]，；（2），x[5,5]. ))1x234x2解 将计算n次拉格朗日插值多项式Ln(x)的值的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: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次拉格朗日插值多项式Ln(x)的图形的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+4x^2)))'),pause n=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 off n=5; x0=-pi:2*pi/(n-1):pi;
y0=tan(cos((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 off n=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 off n=9; x0=-pi:2*pi/(n-1):pi;
86．
第三篇 第六章 函数的插值方法
y0=tan(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 off n=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工作窗口输入名为 Li1.m的M文件的文件名
>> Li1.m
回车运行后，便会逐次画出f(x)在区间[,]上的n次拉格朗日插值多项式Ln(x)
(n2,4,6,8,10)的图形(略).
（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+x^2)'),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=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(x,y5,'m'),gtext('n=10') title('高次拉格朗日插值的振荡')
2回车运行后，便会逐次画出y1/(1x)在区间[5,5]上的n次拉格朗日插值多项式
Ln(x) (n2,4,6,8,10)的图形(略).
6.5.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)','分段线性插值函数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=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);
87．
第三篇 第六章 函数的插值方法
x=-6:0.001:6; y=cos(x); plot(x0,y0,'o',xi,yi,x,y,':'), legend('节点(xi,yi)','分段线性插值函数','被插值函数y=cosx ') title('y=cosx及其分段线性插值函数和节点的图形')
1，在区间[1,1]上取等距节点(xi,yi)i0,1,2,,10, 2125x构造分段线性插值函数Sn(x)，用MATLAB程序计算各小区间的中点xi处Sn(x)的值，
例6.5.5 设函数f(x)作出节点，插值点，f(x)和Sn(x)的图形.
解 节点的横坐标和插值点等取值与例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 x^2)的分段线性插
值的有关图形')
运行后屏幕显示各小区间中点xi处Sn(x)的值，出现节点，插值点，f(x)和Sn(x)的图形(略)
s =
Columns 1 through 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
6.5.5 用MATLAB计算有关分段线性插值的误差
3sin2x)),在区间[,]上取等距节点234x(xi,yi)，i0,1,2,,9，构造分段线性插值函数Sn(x).
(1)用MATLAB程序计算各小区间中点xi处Sn(x)的值，作出节点，插值点，f(x)和
例6.5.8 设函数f(x)tan(cos(Sn(x)的图形；
(2) 用MATLAB程序计算各小区间中点处Sn(x)的值及其相对误差； (3) 用MATLAB程序估计maxf(x)和Sn(x)在区间[,]上的误差限.
x\"解 在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,Ri,i=[xi',fi',si',Ri']
title('y=tan(cos((sqrt(3)+sin(2 x))/(3+4 x^2)))的分段线性插
值的有关图形')
运行后屏幕显示Ri （略），并且作出节点，插值点，f(x)和Sn(x)的图形(略).
在MATLAB工作窗口输入程序
88．
第三篇 第六章 函数的插值方法
>> syms x y=tan(cos((3^(1/2)+sin(2*x))/(3+4*x^2)));yxx=diff(y,x,2),
运行后屏幕显示函数f(x)的二阶导数f''(x)（略）.
在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.*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.*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=(h^2)* myxx/8 运行后屏幕显示myxx =maxf(x)和Sn(x)在区间[,]上的误差限R如下
x\"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 试用MATLAB程序计算例6.6.1中在各小区间中点处分段三次埃尔米特插值Hn(xi1/2)及其相对误差.
解 在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 (略).
6.6.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: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')
title(' 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')
89．
第三篇 第六章 函数的插值方法
title(' y=cosx及其分段埃尔米特插值函数和节点的图形')
例6.6.6 设函数f(x)1定义在区间[5,5]上，节点(X(i),f(X (i)))的横坐1x2标向量X的元素是首项a=-5，末项b=5，公差h=1.5的等差数列，构造三次分段埃尔米特插值函数Hn,3(x).把区间[5,5]分成20等份，构成20个小区间，用MATLAB程序计算各小区间中点xi处Hn,3(x)的值，并作出节点，插值点，f(x)和Hn,3(x)的图形.
解 在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+x^2)及其分段埃尔米特插值函数，插值，节点(xi,yi)
的图形')
运行后屏幕显示各小区间中点xi处Hn,3(x)的值，出现节点，插值点，f(x)和
Hn,3(x)的图形(略).
例6.6.7 设函数f(x)0.5xcosx定义在区间[,]上，取n7，按等距节点构造分段埃尔米特插值函数H7,3(x)，用MATLAB程序计算各小区间中点xi处
H7,3(x)的值，作出节点，插值点，f(x)和H7,3(x)的图形.
解 记节点的横坐标xiih,h2/7,i0,1,2,,7插值点
xi121(xixi1)，i0,1,2,,6.在MATLAB工作窗口输入程序 2>> 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) 的图形')
运行后屏幕显示各小区间中点xi处H7,3(x)的值，出现节点，插值点，f(x)和H7,3(x)的图形(略).
6.6.4 用MATLAB计算有关分段埃尔米特插值的误差
1定义在区间[1,1]上，取n10，按等距节点构2125x(4)造分段埃尔米特插值函数Hn,3(x)，用MATLAB程序在[1,1]上计算maxf(x)和
例6.6.8 设函数f(x)1x1Hn,3(x)的误差公式和误差限.
解 在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=(h^4)* abs(myxxxx/384)
运行后输出f(x)的4阶导数在区间[1,1]上绝对值的最大值myxxxx和Hn,3(x)在区间
[1,1]上的误差公式myxxxx为
myxxxx = R =
15000 625/16*h^4
（4）在MATLAB工作窗口输入程序
90．
第三篇 第六章 函数的插值方法
>> h=0.2; R =625/16*h^4 运行后输出误差限为
R =
0.06250000000000
3sin2x取n9，按等))定义在区间[,]上，234x距节点构造分段埃尔米特插值函数Hn,3(x).
例6.6.9 设函数f(x)tan(cos((1)用MATLAB程序计算各小区间中点xi处Hn,3(x)的值，作出节点，插值点，f(x)和Hn,3(x)的图形；
(2) 并用MATLAB程序计算各小区间中点处Hn,3(x)的值及其相对误差； (3) 用MATLAB程序求maxfx(4)(x)和Hn,3(x)在区间[,]上的误差公式和各
插值的误差限.
解 （1）记节点的横坐标xiih,h2/9,i0,1,2,,9，插值点x1(xixi1)，
i212i0,1,2,,8.在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,fi,Hi,Ri,i=[xi',fi',Hi',Ri']
title('函数y=tan(cos((sqrt(3)+sin(2x))/(3+4x^2)))及其分段埃
尔米特插值函数，插值，节点(xi,yi) 的图形')
运行后屏幕显示各小区间中点xi处的函数值fi，插值Hi，相对误差值Ri，并且作出节点，插值点，f(x)和Hn,3(x)的图形(略).
（2）在MATLAB工作窗口输入程序
>> syms x
y=tan(cos((3^(1/2)+sin(2*x))/(3+4*x^2))); yxxxx=diff(y,x,4),%simple(yxxxx)
运行后屏幕显示函数f(x)的4阶导数f(4)(x)，然后将输出的fx(4)(x)编程求
maxf(4)(x)和Hn,3(x)及其在区间[,]上的误差限的MATLAB程序如下
>>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.*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./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))./(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)).
91．
第三篇 第六章 函数的插值方法
/(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./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(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).*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+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./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))./(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+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((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.*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+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.*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(2.*x)./(3+4
92．
第三篇 第六章 函数的插值方法
.*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+4.*x.^2).^3.*x)
myxx=max(yxxxx), R=(h.^4).* abs(myxx./384)
将其保存为名为myxx.m的M文件，然后在MATLAB工作窗口输入该文件名
>> myxx 运行后屏幕显示myxx =maxfx(4)(x)和Hn,3(x)在区间[,]上的误差公式
h4Rmaxf(4)(x)如下
x384myxx = R =
73.94706841647552 1734520780029061/9007199254740992*h^4
最后在MATLAB工作窗口输入
>>h=2*pi/9; R =1734520780029061/9007199254740992*h^4
运行后屏幕显示Hn,3(x)在区间[,]上的误差限
R =
0.04574453029948
6.7 三次样条（Spline）及其MATLAB程序
6.7.4 用一阶导数计算的几种样条函数和MATLAB程序
计算压紧样条的MATLAB主程序
function
[m,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);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,m)=2; A(i,i+1)=lambda(i); A(i+1,i)=mu(i); end
dY=A\\D'; 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主程序ClampedSp.m求例6.7.1的问题.
93．
第三篇 第六章 函数的插值方法
解 保存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]=ClampedSp(X,Y,dyo,dyn)
运行后输出压紧样条函数sk，X的维数m，由hi1(i1,2,,n)、i、i和di （6.90）式的系数矩阵A，(i1,2,,n1)的坐标组成的向量H、lambda和mu、D，
线性方程组（6.90）的解向量dY如下
sk =
2*(6-2*x)*x^2+2*x*(x-2)^2+9/4*(x-2)*x^2 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-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)*(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
运行后输出f(3)s2(3)和f(8)s4(8)的值如下
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]=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(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
94．
第三篇 第六章 函数的插值方法
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)= 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=A\\D'; 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为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，由hi1(i1,2,,n)、i、i和di
(i1,2,,n1)的坐标组成的向量H、lambda和mu、D，（6.88）式的系数矩阵A，
线性方程组（6.88）的解向量dY如下
sk =
2*(6-2*x)*x^2+915/172*x*(x-2)^2+117/86*(x-2)*x^2 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)*(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 =
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
95．
第三篇 第六章 函数的插值方法
285/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 运行后输出f(3)s2(3)和f(8)s4(8)的值如下
x2 = s2 = x 4= s4 = 3 24.6221 8 68.5581
例6.7.3和例6.7.2用的方法不同，所以计算的结果不同，但是两种方法计算的结果相近.
6.7.6 用MATLAB计算三次样条
例6.7.6 给出节点的数据如下： -2.46 -5.26 -1.87 0.05 1.65 2.69 4.56 7.89 10.31 分别求在下列条件下在插值点x10.02，x22.56处的压紧三次样条插值，并显示该样条函数的有关信息：
(3.50)29.16； (1)5，Sn（1）端点约束条件为Snx -1.00 -0.54 0.13 1.12 1.89 2.06 2.54 2.82 3.50 f(x) (3.50)0. (1)0，Sn（2）端点约束条件为Sn解 （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), PP = spline (X, [5,Y,29.16])
运行后屏幕显示压紧样条分别在x10.02，x22.56处的插值和该样条函数的有关信息如下
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）因为端点约束条件为Sn(1)0，S'n(3.50)0，所以输入MATLAB程序
>> YI= spline (X, [0,Y,0],XI), PP= spline (X, [0,Y,0])
运行后屏幕显示压紧三次样条分别在x10.02，x22.56的插值和该样条函数的有关信息如下
YI =
96．
第三篇 第六章 函数的插值方法
-3.0192 4.7501 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 在默认端点约束条件下，用两种方法计算在下列插值点处的三次样条插值.
（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= 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.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
例6.7.8 设函数f(x)1定义在区间[1,1]上，取n10，按等距节点
125x2构造分段三次样条函数Sn(x).试用MATLAB程序计算在各小区间中点处分段三次样条插值Sn(xi1/2)及其相对误差.
解 .在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(略).
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,'*',x,y,'-.')
legend('节点(xi,yi)', '分段三次样条函数','插值点(x,S)','被插值
函数y')
97．
第三篇 第六章 函数的插值方法
（二）不限定端点约束条件的作图程序
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)', '分段三次样条函数','插值点(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 (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 = 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,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)
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,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))); 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)./y);
RX=abs((y-yX)./y);RLj=abs(y-yL);mRLj=mean(RLj); RSj=abs(y-yS);
98．
第三篇 第六章 函数的插值方法
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'], plot(x0,y0,'bo',X,Y,'k-',X,YS,'rp', X,YH,'g>'),
legend('节点','被插值函数','三次样条函数','分段埃尔米特插值函数')
例6.7.9 设函数f(x)1定义在区间[5,5]上，节点(X(i),f(X (i)))的横坐21x标向量X的元素是首项a=-5，末项b=5，公差h=1.5的等差数列，构造分段三次样条函数Sn(x).把区间[5,5]分成20等份，构成20个小区间，用限定端点约束条件和不限定端点约束条件的MATLAB程序计算各小区间中点xi处Sn(x)的值，并作出节点，插值点，
f(x)和Sn(x)的图形，并与分段埃尔米特插值函数的图形比较.
解 记节点的横坐标xi5ih,h0.5,i0,1,2,,20,插值点xi121(xixi1)，2i0,1,2,,19.
(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+x^2)及其三次样条插值函数,节点和插值点的图形')
运行后屏幕显示各小区间中点xi处Sn(x)的值，作出节点，插值点，f(x)和Sn(x)的图形(略).
'（2）限定端点约束条件，取dy1dynSn(5)0，在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+x^2)及其分段压紧三次样条函数,节点和插值点的图形')
运行后屏幕显示各小区间中点xi处Sn(x)的值（略），作出节点，插值点，f(x)和Sn(x)的图形(略).
如果调节端点约束条件或者增加节点的倍数，例如，在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+x^2)及其增加节点后的不限定端点约束条件的三次样条函
数,节点和插值点的图形')
则运行后输出的图形中的三次样条函数与被插值函数的图像基本重合.
例6.7.10 设函数f(x)0.5xcosx定义在区间[2,2]上，取n7，按等距节点构造分段三次样条函数Sn(x)，用MATLAB程序计算各小区间中点xi处Sn(x)的值，分别作出局部和整体区间上的节点，插值点，f(x)，三次样条函数Sn(x)和分段三次埃尔米特插值函数Hn(x)的图形，并进行比较.
解 编写并保存名为sanci.m的M文件如下
subplot(2,2,1)
99．
第三篇 第六章 函数的插值方法
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');
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, '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,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=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 的函数')
在MATLAB工作窗口输入文件名 >> sanci
运行后屏幕显示各小区间中点xi处三次样条插值pp和分段三次埃尔米特插值Y，出现被插值函数f(x)，节点，三次样条和分段三次埃尔米特的函数及其插值点等图形（略）。
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
6.7.8 用MATLAB计算有关分段三次样条的误差
例6.7.12 设函数f(x)tan(cos(3sin2x))定义在区间[,]上，取n7，234x按等距节点分别作分段线性插值、拉格朗日插值、三次样条插值和分段埃尔米特插值.用MATLAB程序计算各小区间中点xi处四种插值及其绝对误差和相对误差，作出f(x)及其后三种插值函数，插值点和节点的图形，并进行比较.
解 编写并保存名为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;
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)));
100．
第三篇 第六章 函数的插值方法
YL=lagr1(x0,y0,X); YS=interp1(x0,y0,X,'spline'); YH=interp1(x0,y0,X, 'pchip'); yL=lagr1(x0,y0,x); yX=interp1(x0,y0,x); yS=interp1(x0,y0,x,'spline'); yH=interp1(x0,y0,x, 'pchip');
RL=abs((y-yL)./y);RS=abs((y-yS)./y);
RH=abs((y-yH)./y);
RX=abs((y-yX)./y);RLj=abs(y-yL);mRLj=mean(RLj); RSj=abs(y-yS);
mRSj=mean(RSj);RHj=abs(y-yH);RXj=abs(y-yX);mRHj=mean(RHj
);
mRXj=mean(RXj);mRL=mean(RL);mRX=mean(RX); mRS=mean(RS);mRH=mean(RH);
CZ=[x' y' yL' yX' yS' yH'],R=[x' RL' RX' RS' RH'], mR=[mRL' mRX' mRS' 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=tan(cos((sqrt(3)+sin(2x))/(3+4x^2)))及其三种插值
函数,节点和插值点的图形')
运行程序
>> n=7;sanli679
得到各小区间中点x1处的函数值y(x1)，分段线性插值Xn(x1)、拉格朗日插值
i2i12i2i12i2Ln(xi12)、三次样条插值Sn(x)和分段埃尔米特插值Hn(x)及其绝对误差、相对
误差和平均值的结果,作出f(x)及其后三种插值函数，插值点和节点的图形(略).
例6.7.13（机床加工） 待加工零件的外形根据工艺要求由一组数据（x, y）给出（在平面情况下），用程控铣床加工时每一刀只能沿x方向和y方向走非常小的一步，这就需要从已知数据得到加工所要求的步长很小的（x, y）坐标.表 6–15给出的x,y数据位于机翼断面的下轮廓线上（如图6–25），假设需要得到x坐标每改变0.1时的y坐标.试完成加工所需数据，画出曲线，并求出x=0处的曲线斜率和13≤ x ≤15范围内y的最小值.
表 6–15 机翼断面下轮廓线上的部分数据 X Y 0 0 3 1.2 5 1.7 7 2.0 9 2.1 11 2.0 12 1.8 13 1.2 14 1.0 15 1.6 y x Ｏ
图6–25 机翼断面轮廓线（表 6–15数据用圆点表示）
解 根据上述提出的加工要求，以所给数据为节点，在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;
101．
第三篇 第六章 函数的插值方法
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);
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
运行后得到的拉格朗日插值、分段线性插值、三次样条插值和分段埃尔米特插值及其节点的图形，同时还得到拉格朗日插值、分段埃尔米特插值、分段线性插值和三次样条插值的结果(略).
6.8 高元插值及其MATLAB程序
6.8.1 MESHGRID命令的功能和调用格式
调用格式一 [X,Y] =meshgrid (x,y)
例6.8.1 已知x=-3:0.2:3；y=x，计算函数z73xe的值，并作出函数的图形.
解 输入程序
>> [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 X^4 exp(-X^2-Y^2)的图形')
运行后输出函数值和图形（略）.
例6.8.2 作出函数z2xe在区域2x2，2y2上的图形. 解 输入程序
>> [X,Y] = meshgrid(-2:.2:2, -2:.2:2);Z = 2+X .* exp(-X.^2 - Y.^2); mesh(Z)
title('Z = 2+X exp(-X^2 - Y^2)的图形')
运行后输出函数值和图形（略）.
调用格式二 [X,Y] = meshgrid (x)
调用格式三 [X,Y,Z] = meshgrid (x,y,z)
x例6.8.3 计算函数U2xe
24x2y2x2y2y2z2在x0,1,2,4,y0,1,3,5,zy处的函
数值.
解 输入程序
>> x=[0,1,2,-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程序
,-1.5000 , 0 , 1.5000 , 3.0000,yx，例6.8.4 设节点(x,y,z)中的值x -3.00002222x12x2(1y)210(x3y5)exye(1x)y，作在节点(x,y,z)函数z3(1x)e
53 102．
第三篇 第六章 函数的插值方法
处X[2,3,1,7], Y[5,2,-1,5]的双线性插值及其图形.
解 输入程序
>> [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(-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) 的双线性插值图形')
运行后屏幕显示双线性插值及其图形(略).
3xy例6.8.5 设节点(x,y,z)中的x -3:0.5:3,yx和函数Z73xe，作Z在插值点X-3.9:0.5:5, Y-4.9:0.5:4.5处的二元样条插值，双三次插值和数据点的图形.
解 （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 = interp2(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 x^3 exp(-x^2 - y^2) 的二元样条插和值数据点的图形') 运行后屏幕显示Z在插值点X-3.9:0.5:5, Y -4.9:0.5:4.5处的二元样条插值及其图形(略).
（2）计算双三次插值.输入程序
>> [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=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-3x^3exp(-x^2-y^2) 的双三次插值和数据点(x,y,z)的图形')
运行后屏幕显示Z在插值点X-3.9:0.5:5, Y -4.9:0.5:4.5处的双三次插值和数据点图形（略）（三种方法比较留给读者）
4例6.8.6 设节点(x,y,z)中的x -5:0.5:5,yx和函数Z73xe22x2y2，
作Z在插值点X-3.9:0.5:5,Y-4.9:0.5:4.5处的双三次插值和二元最近邻插值及其图形.
解 （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,'r.', 'markersize',3*5) hold off
xlabel('x'), ylabel('y'), zlabel('z'),
title('z =7-3 x^4 exp(-x^2 - y^2) 的双三次插值和数据点的图形')
103．
第三篇 第六章 函数的插值方法
运行后屏幕显示Z在插值点X-3.9:0.5:5,Y -4.9:0.5:4.5处的双三次插值及其图形(略).
（2）二元最近邻插值.输入程序
>> [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, '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 x^3 exp(-x^2-y^2) 的二元最近邻插值和数据点的图形
')
运行后屏幕显示Z在插值点X-3.9:0.5:5,Y -4.9:0.5:4.5处的二元最近邻插值及其图形(略).
6.8.3 三元插值及其MATLAB程序
例6.8.7 设节点(x,y,z)的坐标为x4,0,1,12,y1,0,3,15,zy，计算函数
V2xex2y2z2在插值点xi3:0.25:10,yi3:0.25:3,zi3:0.25:13处的三元线性插值，并作其图形.
解 输入程序
>> 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,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(-x^2 - y^2-z^2) 的三元线性插值图形') hold on,colorbar('horiz'), view([-30 45])
运行后屏幕显示三元线性插值及其图形（略）.
例6.8.8 取 n=10 ，作函数flow在插值点xi0.1:0.25:10,yi3:0.25:3,
ziyi处的三元三次样条插值及其图形.
解 输入程序
>> [x,y,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')
运行后屏幕显示三元三次样条插值及其图形（略）.
104．
因篇幅问题不能全部显示，请点此查看更多更全内容