2013-2014(1)专业课程实践论文
题目:埃尔米特插值
一、算法理论
1、埃尔米特插值多项式
设已知函数yf(x)在节点x0x1xn上的函数值要求一个插值多yif(xi)(i0,1,,n)以及一切导数值yif(xi)(i0,1,n),项式H(x),使其满足
H(xi)yi,H(xi)yi ,i0,1,,n (1) 显然,由条件(1)可以确定一个次数不高于2n1的代数多项式H2n1(x),曲线yH2n1(x)与yf(x)在节点处不仅重合而且有公共切线。我们采用拉格朗日插值基函数的方法。先求插值基函数j(x),j(x)(j0,1,,n),共2n2个基函数,每一个基函数都是一个2n1次多项式,且满足条件
j(x)jk,j(xk)0;j(xk)0,j(xk)jk,j0,1,,n.
这里
0,jk,
1,jk.
n
(2)
(3)
于是满足条件(1)的插值多项式 可写成用插值函数表示的形式
H2n1(x)[yj(x)yjj(x)]
j0
(4)
n1(xi)yj(i0,1,,n).下面的问题就由条件(2),显然有H2n1(xi)yi,H2
是要求满足条件(2)的j(x)与j(x).为此,可利用拉格朗日插值基函数lj(x),由条件(2)有n个二重零点xk(k0,1,,n,kj),于是可令
j(x)(axb)l2j(x). 由条件(2)有
j(xj)(axb)l2j(xj)1,
j(xj)lj(xj)[alj(xj)2(axjb)lj(xj)]0.解出
a2lj(xj),b12xjlj(xj). 由于
lj(x)
xxj1xxj1xx0xx1xxn
,
xjx0xjx1xjxj1xjxj1xjxn
n
故
lj(xj)
k0kj
1
, xjxk
n
于是
j(x)[12(xxj)
k0kj
1
]l2j(x).xjxk
(5)
同理,可得
j(x)(xxj)l2j(x). (6) 将(5)式、(6)式 代入式(4)便得到埃尔米特插值多项式
n
12
H2n1(x)yi[12(xxi)]lj(x)yj(xxj)l2j(x)
xxj0k0jj0k
kj
n
n
(7)
满足条件(1)的埃尔米特插值多项式是唯一的。这可用反证法证明,此处
从略。
2、两点三次埃尔米特插值
设已知yf(x)在[a,b]上的节点x0,x1上的函数值y0,y1及一阶导数值
,y1,则可按公式(7)写出三次埃尔米特插值多项式 y0
0(x)y11(x)H3(x)y00(x)y11(x)y0 y0(12
xx0xx12xx1xx02
)()y1(12)() x0x1x0x1x1x0x1x0
1
y(xx)(xx
x0x1
(xx1)()2y1
xx02
).x1x0
二、算法框图
两点三次埃尔米特插值框图如下:
两点三次埃尔米特插值框图如下:
1.实现埃尔米特插值的MATLAB函数文件hermite2.m如下
function yy=hermite2(x,y,dy,xx)
n=length(y);m=length(x);l=length(dy);k=length(xx);
if m~=n,error('向量长度不一致');end;
if n~=l,error('向量长度不一致');end;
z=zeros(1,k);
for j=1:k
s=0;
for t=1:m;
a=0;b=1;
for i=1:n;
if x(t)~=x(i)
a=a+1/(x(t)-x(i));
b=b*((xx(j)-x(i))/(x(t)-x(i)));
end
end
s=s+(y(t)*(1-2*(xx(j)-x(t))*a)*b^2+dy(t)*(xx(j)-x(t))*b^2); end
z(j)=s;
end
yy=z;
2.实现两点三次埃尔米特插值的MATLAB函数文件hermite.m如下 function yy=hermite(x,y,dy,xx)
n=length(y);m=length(x);l=length(dy);k=length(xx);
if m~=n,error('向量长度不一致');end;
if n~=l,error('向量长度不一致');end;
z=zeros(1,k);
for i=1:k;
s=0;
a1=(1-2*(xx(i)-x(1)/(x(1)-x(2))))*((xx(i)-x(2))/(x(1)-x(2)))^2; a2=(1-2*(xx(i)-x(2)/(x(2)-x(1))))*((xx(i)-x(1))/(x(2)-x(1)))^2; b1=(xx(i)-x(1))*((xx(i)-x(2))/(x(1)-x(2)))^2;
b2=(xx(i)-x(2))*((xx(i)-x(1))/(x(2)-x(1)))^2;
s=y(1)*a1+y(2)*a2+dy(1)*b1+dy(2)*b2;
z(i)=s;
end
yy=z;
例1.设次数不超过3的多项式P3(x)满足插值条件
P,P,P3(0)0,P3(0)13(1)13(1)2
求其在x=0.5,0.75处的近似值。
解:在MATLAB命令窗口输入:
format long;x=[0;1];y=[0;1];dy=[1;2];xx=[0.5;0.75];
yy=hermite(x,y,dy,xx)
其得到结果如下:
yy =
0.[**************] 0.[**************]
例2.设次数不超过5的多项式P5(x)满足插值条件
P,P,P3(0)0,P3(0)13(1)13(1)2,P3(2)3,P3(2)5
求其在x=0.5,0.75,1.15处的近似值。
解:在MATLAB命令窗口输入:
format long;x=[0;1;2];y=[0;1;3];dy=[1;2;5];xx=[0.5;0.75;1.15]; yy=hermite2(x,y,dy,xx)
其得到结果如下:
yy =
0.[**************] 0.[**************] 1.[**************]