feature-lima
6.1
插值问题及其误差
6.1.2
与插值有关的
MATLAB
函数
(
一
)
POLY2SYM
函数
调用格式一
：
poly2sym
(C)
调用格式二
：
f1=poly
2sym(C,'V')
或
f2=poly2sym(C, sym ('V')
)
,
(
二
)
POLYVAL
函数
调用格式
：
Y =
polyval(P
,X)
(
三
)
POLY
函数
调用格式
：
Y = poly
(V)
(
四
)
CONV
函数
调用格式
：
C =conv (A,
B)
例
6.1.2
求三个一次多项式
、<
/p>
和
的积
.
它们的
零点
分别依次为
0.4,0.8,1.2.
解
我们可
以用两种
MATLAB
程序求之
.
方法
1
<
/p>
如输入
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)
(
六
)
ro
ots(poly(1:n))
命令
调用格式
：
roots(poly(1:n))
(
七
)
det(a*eye(size (A)) -
A)
命令
调用格式
< br>：
b=det(a*ey
e(size (A)) -
A)
6.2
拉格朗日（
Lagrange
）插值及其
MATLAB<
/p>
程序
6.2.1
线性插值及其
MATL
AB
程序
例
6.2.1
已知函数
在
上具有二阶连续导数，
，并估计其误差
.
，且满足条件
.<
/p>
求线性插值多项式和函数值
解
输入程序
>>
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
求函数
e
在
上线性插值多项式，并估计其误差
.
解
输入程序
>> 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 =
-14234/22548*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]
< br>分成
等份
，
用
< br>产生
个
.
节点，然后根据（
6.9
）和（
6.13
）式分别作线性插值函数
和抛物线插值函数
用它们分别计算
cos (
π
/6) (
取四位有效数字
)
，并估计其误差
< br>.
解
输入程序
>>
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
< br>、插值多项式和插值为
l0 =
-57349/9992*x+1
l1 =
57349/9992*x
P =
-0.6366 1.0000
L =
-57349/9992*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 =
2285/2856*x^2-2157/42624*x+1
l1 =
-2285/148*x^2+57349/22548*x
l2 =
2285/2856*x^2-57349/9992*x
P
=
-0.3357 -0.1092 1.0000
L=
-6875/18984*x^2-7879/7293
6*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
例
6.2.4
给
出
节
点
数
据
，
，
，
次拉格朗日（
La
grange
）插值及其
MATLAB
程序
，作三次拉格朗日插值多项式计算
，并估计其误差
.
解
输入程序
>>
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),
运行后输出基函数<
/p>
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
，
.
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
给出节点数据
，
，
，
，
，
，作五次拉格朗日插值多项式和基函数，
并写出估计其误差的公式
.
解
在
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<
/p>
如下
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+<
/p>
0.0648
*x^4
-0.2169<
/p>
*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.050
3*x^3+0.6305*x^2-0.4852*x+0.0048]
[
-0.0693*x^5+0.2
184*x^4+0.3961*x^3-1.2116*x^2-0.3166*x+1.0033]
[
<
/p>
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.25
30*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]
估计其误差的公式为
，
.
6.2.6
拉格朗日插值及其误差估计的
< br>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
已知
，<
/p>
，
，用拉格朗日
插值及其误差估计的
MATLAB
主程序求
的近似值，并估计其
误差
.
解
在
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)
运行后输出插值
< br>y
及其误差限
R
为
y = R =
0.6434 8.8610e-004.
6.3
牛顿（
< br>Newton
）插值及其
MATLAB
< br>程序
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
给
出
节
点
数
据
，
，
，
，
，
,
作五阶牛顿插值多项式和差商，并写出
其估计误差
的公式
.
解
（
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
及其系数向量<
/p>
C
,
插值余项公式
L
及其向量
C
w
< br>如下
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.1
129
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 =
-78629/36968*x^5+583849564517807/900719925
4740992*x^4+593245028711263/2856*x^3+3823
593773002357/11258999
06842624*x^2-3219/
74*x+399/28147497671
0656
wcgs =
1/720*M*(x^6-79/25*x^
5-14201/2500*x^4+4934/28
56*x^3+1545/352
*x^2-8179/56294995342131
2*x+521276/3696
8)
Cw =
0.0014
-0.0044
-0.0079
0.0243
0.0061
-0.0202
0.0002
即
L =1.0954-4.574
5*
x
+3.3960*
x
^2+2.1076*
x
^3+0.0648*<
/p>
x
^4-0.2169*
x
^5.
估计其误差的公式为
，
例
6.3.4
求函数
e
在
上五阶牛顿插值多项式
，估计其误差的公式
.
和误差限公式
.
用它们计算
，并估计其误差
.
解
（
1
）输入
MATLAB
程序
< br>
>> 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/11529215*x
^5-3515/922337
2*x^4+169/18984*x^3-53501
/
9992*x^2+62987/4596*x-394
3/
562949953421312
wcgs =
1/720
*M*(x^6-24*x^5+1172/5*x^4-5952/5*x^3+7276634802928
539/219902
3255552*x^2-5237461
186650519/1*x+6121/21990232555
52)
Cw =
0.0014
-0.0333
0.3256
-1.6533
4.5959
-6.6159
3.8438
M =
-1.3494e-004
（
2
）输入
MATLAB
程序
>> syms x
wc
gs1=1/720*M*1/720*M*(x^6-24*x^5+1172/5*x^4-5952/5*
x^3+7
276634802928539/22*x^2-52374611866
50519/1*
x+6121/22),
运行后输出误差限
公式
wcgs1
如下
wcgs1 =
5565367633581443/790
0672*x^6-16696102
900744329/987584*x^5+1
632799/1
98843983859875840*x^4-5074199/8
53840*x^3+45501777/3484491437
27598944*x^2-29471
8917/3529472*x+3387<
/p>
83926665276375603/3484494
（
3
）输入
MATL
AB
程序
>> x=3.156;
y=9721799720875/11529215*x^5-3515/9
223372*x^4+169/18984*x^3
-1251
1/9992*x^2+62987/45
96*x-3943/56
2949953421312
wcgs2=1/720*M*(x^6-24*x^5
+1172/5*x^4-5952/5*x^3+7276634802
928539
/22*x^2-5237461186650519/1*x+6085939
356
447121/22)
运行后输出
的近似值
< br>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
已知
求
，
，
，用牛顿插值法
的近似值，估计其误差，并与例
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 =
-08357/4596*x^2+
1405/11
25899906842624*x-1324/22548
wcgs =
1/6*M*(x^3-3/4*x^2*pi
+4799/22548*x-
7757769783530263/18984)
Cw =
0.1667 -0.3927
0.2970 -0.0718
y =
0.6434
上述两种方法计算
y
的
结果相同
.
6.3.5
牛顿插值法的
MATL
AB
综合程序
求牛顿插值多项式、差
商、插值及其误差估计的
MATLAB
主程序
< br>
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
给出节点数据
，
，
，
，作三阶牛顿插值多项式，计算
，并估计其误差
.
解
首先将名为
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
及其误差限公式
R
，
三阶牛顿插值多项式
< br>P
及其系数向量
C
，
差商的矩阵
A
如下
y =
22.3211
R =
1323/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]
分成
等份
，用
产生
个节
点，求二阶和三阶牛顿插值多项式
和
.
用它们分别计算
sin
(
π
/7) (
取四位有效数
字
)
，并估计其误差<
/p>
.
解
首先将名为
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
=
和
y
2
=<
/p>
及其误差限
R
1
和
R
2
，
二阶
和三阶牛顿插值多项式
P
2
和
P
3
及其系数向量
C<
/p>
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 =
-3437/9992*x^2+16382/148*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 =
-1963/9992*x^3-47127/72
0
57594037927936*x^2+45956/4596*x
6.4
埃尔米特（<
/p>
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
给
定
函
数
，
在
点
，
处
的
函
数
值
，
,
，
和误差公
和导数值
在点
且
，求函数
处的五阶埃尔米特插值多项式<
/p>
式，计算
并估计其误差
.
解
（
1
）保存名为
hermite.m
的
M
文件
.
（
2
）在
MATLAB<
/p>
工作窗口输入程序
>>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/352*x^5-4209/439804
6511104*x^4+7434/4398046511104*x^3-31622/2
2*x^2+5835/8796093022208*x-684
3/74
wcgs =
1/720*M*(x^6-11/6*x^5
*pi+74467/562949953421312*x
^4-436374550
3235773/2856*x^3+255/2199023255
552*x^2-
7178/22548*x+375843/9007
2)
Cw =
0.0014 -0.0080
0.0184 -0.0215 0.0136 -0.0044 0.0006
当
时的误差公式为
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<
/p>
工作窗口输入程序
>>
x=1.567;M=1;
Hk=6725828781679091/352*x
^5-4209/43
98046511104*x^4+7434/43980465
11104*x^3-31622
3/22*x^2+5835/8796093022
208*x-6
843/74,
wcgs=1/720*M*
(x^6-11/6*x^5*pi+74467/56294995342
1312*
x^4-4363745503235773/2856*x^3+255/21990
23255552*x^2-7178/22548*x+375843/
9992)
运行后输出
的近似值
Hk
及其误差
wcgs
如下
Hk =
2.5265
wcgs =
1.3313e-008
6.5
分段插值及其
MATLAB
程序
6.5.1
高次插值的振荡和
MA
TLAB
程序
例
6.5.1
作
下
列
函
数
在
指
定
区
间
上
的
的图形，并讨论插值多项式
次
拉
格
朗
日
插
值
多
项
式<
/p>
的次数与误差
的关系
.
（
1
）
，
；
（
2
）
，
.
解
将计算
次拉格朗日插值多项式
的值的
MATL
AB
程序保存名为
lagr1.m
的<
/p>
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
）
现提供作
次拉格朗
日插值多项式
的图形的
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-'
), <
/p>
gtext(
'y=tan(cos((sqrt(3)+sin
(2x))/(3+4x^2)))'
),pause
n=3; x0=-pi:2*pi/(3-1):pi;
y
0=tan(cos((3^(1/2)+sin(2*x0))./(3+4*x0.^2)));
y1=lagr1(x0,y0,x);hold on,
p
lot(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,
p
lot(x,y2,
'b:'
),gtext(
'n=4'
),pause,hold off
n=7; x0=-pi:2*pi/(n-1):pi;
y
0=tan(cos((3^(1/2)+sin(2*x0))./(3+4*x0.^2)));
y3=lagr1(x0,y0,x);hold on,
p
lot(x,y3,
'rp'
),gtext(
'n=6'
),pause,hold off
n=9; x0=-pi:2*pi/(n-1):pi;
y
0=tan(cos((3^(1/2)+sin(2*x0))./(3+4*x0.^2)));
y4=lagr1(x0,y0,x);hold on,
p
lot(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(
'
高次拉格朗日插值的振荡
< br>'
)
在
MA
< br>TLAB
工作窗口输入名为
L
i1.m
的
M
文件的文件名
>> Li1.m
回车运行后，便会逐次画
出
的图形
(<
/p>
略
)
.
（
2
）
在
MA
TLAB
工作窗口输入程序
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'
)
t
itle(
'
高次拉格朗日插值的振荡
'
)
回车运行后，便会逐次画出
<
/p>
的图形
(
略
)<
/p>
.
在区间<
/p>
上的
次拉格朗日插值多项式
在区间
上的
次拉格朗日插值多项
式
< br>
6.5.4
作有关分段
线性插值图形的
MATLAB
程序
<
/p>
作有关分段线性插值图形的
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,
'-'
,x
i,s,
'*'
,x,y,
'-.'<
/p>
)
legend(
'
< br>节点
(xi,yi)'
,
'
分段线性插值函数
Sn(x)'
,
'
插值点
(x,s)'
,
'
被
插值函数
y'
)
我们也可以直接在在
MA<
/p>
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,
':'
),
lege
nd(
'
节点
(xi,yi)'
,
'
分段线性插值函数
'
,
'
被插值函数
< br>y=cosx '
)
title(
'y=cosx
及其分段线性插值函数和节点的图形
')<
/p>
例
6.5.5
设函数
构造分段线性插值函数
，在区间
上取等距节点
处
,
的值，作
，用
MA
TLAB
程序计算各小区间的中点
出节点，插值点
，
和
的图形
.
解
节点的横坐标和插值点等取值与例
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)
的分段线性插值
的有关图形
')
运行后屏幕显示各小区间中点
处
的值，出现节点，插值点，
和
的图形
(
略
)
s =
Columns 1 through 4
0.65 0.88 0.15
0.35
Columns 5 through 8
0.75
0.75 0.35
0.15
Columns 9 through 10
0.88
0.65
6.5.5
用
MATLAB
计算有关分段线性插值的误差<
/p>
例
6.5.8
设函数
，
，
构造分段线性插值函数
处
,
在区间
.
上取等距节点
(1)
< br>用
MA
TLAB
程序计算各小区
间中点
的值，作出节点，插值点，
和
的
图形；
(2)
用
< br>MA
TLAB
程序计算各小区间中点处
< br>的值及其相对误差；
(3)
用
MA
TLAB
程序估计
和
在区间
上的误差限
.
解
在
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)))
的分段线性插值
的有关图形<
/p>
'
)
运行后屏幕显示
< br>R
i
（略）
< br>，并且作出节点，插值点，
和
的图形
(
略
).
在
MATLAB
工作窗口输入程序
>> syms x
y=tan(cos((3^(1
/2)+sin(2*x))/(3+4*x^2)));yxx=diff(y,x,2),
运行后屏幕显示函数
的二阶导数
（略）
.
在
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(co
s
((3.^(1/2)+sin(2.*x))./(3+4.*x.^2))).^2).*sin((3^(
1/2)+sin(2.*x)
)./(3+4.*x.^2)).^2.*(2*co
s(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)./
< br>(3+4.*x.^2)-8*(3^(1/2)+sin(2.*x))./(3+4.*x.^2). ^2.*x).^2-(1+tan(c
os((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
=
和
在区间
上的误差限
如
下
myxx =
R =
-0.64 -0.11
6.6
分段埃尔米特插值及其<
/p>
MATLAB
程序
6.6.2
分段埃尔米特插值的
M
ATLAB
程序
调用格式一：
YI=interp1(X,Y,XI,'pchip')
调用格式
二
：
YI=interp1(X,XI,'pchip')
例
6.6.5
试用
MA
T
LAB
程序计算例
6.6.1
中在各小
区间中点处分段三次埃尔米特插
值
及其相对误差
.
解
在
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<
/p>
(
略
).
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,
'*'<
/p>
,x,y,
'-.'
)
legend(
'
节点
(xi
,yi)'
,
'
分段埃尔米特插值函
数
'
,
'
插值
点
(x,H)'
,
'
< br>被插
值函数
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)'
,
'
分段埃尔米特插值函数
'
,
< br>'
被插值函数
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
设函数
定义在区间
< br>上，节点
(
X
(
i
),
(
X
< br> (
i
)))
的横坐标
向量
X
的元素是首项
a
=-5
，
末项
b
=5
，
公差
h
=1.5
的等差数列，
构造三次分段
埃尔米特插值
函数
.
把区间
分成
20
等份，构成
20
个小区间，用
MA
TLAB
程序计算各小
区间中点
处
的值，并作出节点，插值点，
和
的图形
.
解
<
/p>
在
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)
的图形
'
)
运行后屏幕显示各小区间中点
的图形
(
略
).
例
6.6.7
设函数
构造分段埃尔米特插值函数<
/p>
值，作出节点，插值点，
和
处
的值，
出现节点，
插值点，
和
定义在区间
上，取
，按等距节点<
/p>
处
的
，用
MA<
/p>
TLAB
程序计算各小区间中点
的图形<
/p>
.
解
记
节
点
的
横
坐
标
插
值
点
，
.
在
MA
TLAB
工作窗口输入程序
>> 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)
的图形
'
)
feature-lima
feature-lima
feature-lima
feature-lima
feature-lima
feature-lima
feature-lima
feature-lima
-
上一篇：刀片材质分类
下一篇：20秋学期《大学英语(四)》平时作业1