第六章函数的插值方法 (2)

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 求三个一次多项式)(x g 、)(x h 和)(x f 的积)()()(x h x g x f ⋅⋅.它们的零点分别依次为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 已知函数)(x f 在]3,1[上具有二阶连续导数，5)("≤x f ，且满足条件2)3(,1)1(==f f .求线性插值多项式和函数值)5.1(f ，并估计其误差.
解 输入程序
>> 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)
运行后输出基函数l 0和l 1及其插值多项式的系数向量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 =
1.8750
例6.2.2 求函数=)(x f e x -在]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),
运行后输出基函数l 0和l 1及其插值多项式的系数向量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.
6.2.2 抛物线插值及其MATLAB 程序
例 6.2.3 求将区间 [0, π/2] 分成n 等份)2,1(=n ，用x x f y cos )(==产生1+n 个节点，然后根据（6.9）和（6.13）式分别作线性插值函数)(1x P 和抛物线插值函数)(2x P .用它们分别计算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)
运行后输出基函数l 0和l 1及其插值多项式的系数向量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) ,
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)
运行后输出基函数l 01、l 11和l 21及其插值多项式的系数向量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 给出节点数据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];
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),
运行后输出基函数l 0，l 1，l 2和l 3及其插值多项式的系数向量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.
输入程序
>> 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
即
R 3 )(04.0)4(ξ≈f ， )00.2,00.2(-∈ξ.
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 给出节点数据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 ，作五次拉格朗日插值多项式和基函数，并写出估计其误差的公式.
解 在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 及其系数矩阵L 1如下
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]
估计其误差的公式为
)(5x R )25.3)(03.2)(02.1)(01.0()00.1)(15.2(!
6)()6(----++=x x x x x x f ξ，)3.25,-2.15(∈ξ.
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 已知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]=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 给出节点数据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 ,作五阶牛顿插值多项式和差商，并写出其估计误差的公式.
解 （1）保存名为newpoly.m 的M 文件.
（2）输入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];
[A,C,L,wcgs,Cw]= newdcw (X,Y)
运行后输出差商矩阵A ,五阶牛顿插值多项式L 及其系数向量C , 插值余项公式L 及其向量C w 如下
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.
估计其误差的公式为
)(5x R )25.3)(03.2)(02.1)(01.0()00.1)(15.2(!
6)()6(----++=x x x x x x f ξ，)3.25,-2.15(∈ξ.
例6.3.4 求函数7)(-=x f e 5/x -在]6,2[上五阶牛顿插值多项式，估计其误差的公式和误差限公式.用它们计算)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)/(5^6)),
运行后输出差商矩阵A , 五阶牛顿插值多项式L 及其系数向量C , 插值余项公式L 及其向量C w 如下
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
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-166961
029********/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)
运行后输出)156.3(f 的近似值y ，及其误差限wcgs 2如下
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 已知5.030sin = ，1707.045sin = ，0866.060sin =
，用牛顿插值
法求
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,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 给出节点数据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,C,P]=newdscg(X,Y,x,M)
运行后输出插值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.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 等份)3,2(=n ，用x x f y sin )(==产生1+n 个节点，求二阶和三阶牛顿插值多项式)(2x P 和)(3x P .用它们分别计算sin (π/7) (取四位有效数字)，并估计其误差.
解 首先将名为newdscg.m 的程序保存为M 文件，然后在MA TLAB 工作窗口输入程序
>> 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)
运行后输出插值y 1=)7/sin()7/(2π≈πP 和y 2=)7/sin()7/(3π≈πP 及其误差限R 1和R 2，二阶和三阶牛顿插值多项式P 2和P 3及其系数向量C 1和C 2，差商的矩阵A 1和A 2如下
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 程序
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 给定函数)(x f 在点2/,4/,6/210π=π=π=x x x 处的函数值5.0)(0=x f , 1707.0)(1=x f ，1)(2=x f 和导数值0866.0)(0'=x f ，1707.0)(1'=x f ，
0000.0)(0'=x f ，且1)()6(≤x f ，求函数)(x f 在点210,,x x x 处的五阶埃尔米特插值
多项式)(5x H 和误差公式，计算)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)
运行后输出五阶埃尔米特插值多项式H k 及其系数向量H c ，误差公式wcgs 及其系数向量C w 如下
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 当M x f =≤1)()6(时的误差公式为
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
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)
运行后输出)1.567(f 的近似值Hk 及其误差wcgs 如下
Hk =
2.5265
wcgs =
1.3313e-008
6.5 分段插值及其MATLAB 程序
6.5.1 高次插值的振荡和MATLAB 程序
例6.5.1 作下列函数在指定区间上的n 次拉格朗日插值多项式)(x L n
)10,8,6,4,2(=n 的图形，并讨论插值多项式)(x L n 的次数与误差)(x R n 的关系.
（1）))432sin 3tan(cos()(2x x x f ++=，],[ππ-∈x ；（2）211)(x
x f +=，]5,5[-∈x . 解 将计算n 次拉格朗日插值多项式)(x L n 的值的MATLAB 程序保存名为
lagr1.m 的M 文件.
function y=lagr1(x0,y0,x)
n=length(x0); m=length(x);
for i=1:m
z=x(i);s=0.0;
for k=1:n
p=1.0;
for j=1: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 次拉格朗日插值多项式)(x L n 的图形的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;
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('高次拉格朗日插值的振荡')
在MA TLAB 工作窗口输入名为 Li1.m 的M 文件的文件名
>> Li1.m
回车运行后，便会逐次画出)(x f 在区间],[ππ-上的n 次拉格朗日插值多项式)(x L n )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+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('高次拉格朗日插值的振荡')
回车运行后，便会逐次画出)1/(12x y +=在区间]5,5[-上的n 次拉格朗日插值多项式
)(x L n )10,8,6,4,2(=n 的图形(略).
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')
我们也可以直接在在MA TLAB 工作窗口编程序.例如，
>> 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);
x=-6:0.001:6; y=cos(x); plot(x0,y0,'o',xi,yi,x,y,':'),
legend('节点(xi,yi)','分段线性插值函数','被插值函数y=cosx ')
title('y=cosx 及其分段线性插值函数和节点的图形')
例6.5.5 设函数2
2511)(x x f +=，在区间]1,1[-上取等距节点),(i i y x 10,,2,1,0 =i , 构造分段线性插值函数)(x S n ，用MA TLAB 程序计算各小区间的中点i x 处)(x S n 的值，作出节点，插值点，)(x f 和)(x S n 的图形.
解 节点的横坐标和插值点等取值与例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)的分段线性插
值的有关图形')
运行后屏幕显示各小区间中点i x 处)(x S n 的值，出现节点，插值点，)(x f 和)(x S n 的图形(略)
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 计算有关分段线性插值的误差
例6.5.8 设函数))432sin 3tan(cos()(2
x x x f ++=,在区间],[ππ-上取等距节点),(i i y x ，9,,2,1,0 =i ，构造分段线性插值函数)(x S n .
(1)用MA TLAB 程序计算各小区间中点i x 处)(x S n 的值，作出节点，插值点，)(x f 和)(x S n 的图形；
(2) 用MA TLAB 程序计算各小区间中点处)(x S n 的值及其相对误差；
(3) 用MA TLAB 程序估计)(max "
x f x ππ≤≤-和)(x S n 在区间],[ππ-上的误差限. 解 在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)))的分段线性插
值的有关图形')
运行后屏幕显示R i （略），并且作出节点，插值点，)(x f 和)(x S n 的图形(略).
在MATLAB 工作窗口输入程序
>> syms x
y=tan(cos((3^(1/2)+sin(2*x))/(3+4*x^2)));yxx=diff(y,x,2),
运行后屏幕显示函数)(x f 的二阶导数)('
'x f （略）.
在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.*c os(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 =)(max "x f x π≤≤π-和)(x S n 在区间],[ππ-上的误差限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/1+i n x H 及其相对误差.
解 在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']
运行后屏幕显示各小区间中点x i 处的函数值f i ，插值s i ，相对误差值R i (略).
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')
title(' y=cosx 及其分段埃尔米特插值函数和节点的图形')
例6.6.6 设函数211)(x
x f +=定义在区间]5,5[-上，节点(X (i ),f (X (i )))的横坐标向量X 的元素是首项a =-5，末项b =5，公差h =1.5的等差数列，构造三次分段埃尔米特插值函数)(3,x H n .把区间]5,5[-分成20等份，构成20个小区间，用MA TLAB 程序计算各小区间中点i x 处)(3,x H n 的值，并作出节点，插值点，)(x f 和)(3,x H n 的图形.
解 在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)
的图形')
运行后屏幕显示各小区间中点i x 处)(3,x H n 的值，出现节点，插值点，)(x f 和
)(3,x H n 的图形(略).
例6.6.7 设函数x x x f cos 5.0)(-=定义在区间],[ππ-上，取7=n ，按等距节点构造分段埃尔米特插值函数)(3,7x H ，用MATLAB 程序计算各小区间中点i x 处)(3,7x H 的值，作出节点，插值点，)(x f 和)(3,7x H 的图形.
解 记节点的横坐标7,,2,1,0,7/2, =π=+π-=i h ih x i 插值点)(2
1121+++=i i i x x x ，6,,2,1,0 =i .在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(x)及其分段埃尔米特插值函数，插值，节点
(xi,yi) 的图形')
运行后屏幕显示各小区间中点i x 处)(3,7x H 的值，出现节点，插值点，)(x f 和)(3,7x H 的图形(略).
6.6.4 用MATLAB 计算有关分段埃尔米特插值的误差
例6.6.8 设函数22511)(x
x f +=定义在区间]1,1[-上，取10=n ，按等距节点构造分段埃尔米特插值函数)(3,x H n ，用MATLAB 程序在]1,1[-上计算)(max )4(1
1x f x ≤≤-和)(3,x H n 的误差公式和误差限.
解 在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)
运行后输出)(x f 的4阶导数在区间]1,1[-上绝对值的最大值myxxxx 和)(3,x H n 在区间
]1,1[-上的误差公式myxxxx 为
myxxxx = R =
15000 625/16*h^4
（4）在MATLAB 工作窗口输入程序
>> h=0.2; R =625/16*h^4
运行后输出误差限为
R =
0.06250000000000
例6.6.9 设函数))432sin 3tan(cos(
)(2
x x x f ++=定义在区间],[ππ-上，取9=n ，按等距节点构造分段埃尔米特插值函数)(3,x H n . (1)用MA TLAB 程序计算各小区间中点i x 处)(3,x H n 的值，作出节点，插值点，)(x f 和)(3,x H n 的图形；
(2) 并用MATLAB 程序计算各小区间中点处)(3,x H n 的值及其相对误差；
(3) 用MA TLAB 程序求)(max )4(x f x ππ≤≤-和)(3,x H n 在区间],[ππ-上的误差公式和各
插值的误差限. 解 （1）记节点的横坐标9,,2,1,0,9/2, =π=+π-=i h ih x i ，插值点)(21
12
1+++=i i i x x x ，8,,2,1,0 =i .在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= hermite tx (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) 的图形')
运行后屏幕显示各小区间中点x i 处的函数值f i ，插值H i ，相对误差值R i ，并且作出节点，插值点，)(x f 和)(3,x H n 的图形(略).
（2）在MATLAB 工作窗口输入程序
>> syms x
y=tan(cos((3^(1/2)+sin(2*x))/(3+4*x^2)));
yxxxx=diff(y,x,4),%simple(yxxxx)
运行后屏幕显示函数)(x f 的4阶导数)()4(x f ，然后将输出的)()4(x f
编程求)(max )4(x f x ππ≤≤-和)(3,x H n 及其在区间],[ππ-上的误差限的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.*c os(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.*t an(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。