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 求三个一次多项式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)
2.2 拉格朗日(Lagrange)插值及其MATLAB程序
2.2.1 线性插值及其MATLAB程序
例2.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为
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 =
76.
1.8750
例2.2.2 求函数f(x)?e?x在[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.
2.2.2 抛物线插值及其MATLAB程序
例2.2.3 求将区间 [0, π/2] 分成n等份(n?1,2),用y?f(x)?cosx产生
n?1个节点,然后根据(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) , 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),
77.
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.
2.2.3 n次拉格朗日(Lagrange)插值及其MATLAB程序
例2.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.
输入程序
>> 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
即
78.
R3 ?0.04f(4)(?), ??(?2.00,2.00).
2.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
例2.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)?(x?2.15)(x?1.00)(x?0.01)(x?1.02)(x?2.03)(x?3.25),??(-2.15,3.25).
6!
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
79.
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 已知sin30?0.5,sin45?0.7071,sin60?0.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.
?2.53 分段插值及其MATLAB程序
2.3.1 高次插值的振荡和MATLAB程序
例2.3.1 作下列函数在指定区间上的n次拉格朗日插值多项式Ln(x) (n?2,4,6,8,10)的图形,并讨论插值多项式Ln(x)的次数与误差Rn(x)的关系.
(1)f(x)?tan(cos(3?sin2x1)),;(2),x?[?5,5]. x?[??,?]f(x)?223?4x1?x解 将计算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;
80.