埃尔米特插值

2013-2014（1）专业课程实践论文

题目：埃尔米特插值

一、算法理论

1、埃尔米特插值多项式

设已知函数yf(x)在节点x0x1xn上的函数值要求一个插值多yif(xi)(i0,1,,n)以及一切导数值yif(xi)(i0,1,n)，项式H(x)，使其满足

H(xi)yi，H(xi)yi ，i0,1,,n (1) 显然，由条件(1)可以确定一个次数不高于2n1的代数多项式H2n1(x)，曲线yH2n1(x)与yf(x)在节点处不仅重合而且有公共切线。我们采用拉格朗日插值基函数的方法。先求插值基函数j(x),j(x)(j0,1,,n)，共2n2个基函数，每一个基函数都是一个2n1次多项式，且满足条件

j(x)jk,j(xk)0;j(xk)0,j(xk)jk,j0,1,,n.

这里

0,jk,

1,jk.

n

(2)



(3)

于是满足条件(1)的插值多项式 可写成用插值函数表示的形式

H2n1(x)[yj(x)yjj(x)]

j0

(4)

n1(xi)yj(i0,1,,n).下面的问题就由条件(2)，显然有H2n1(xi)yi,H2

是要求满足条件(2)的j(x)与j(x).为此，可利用拉格朗日插值基函数lj(x)，由条件(2)有n个二重零点xk(k0,1,,n,kj)，于是可令

j(x)(axb)l2j(x). 由条件(2)有

j(xj)(axb)l2j(xj)1,

j(xj)lj(xj)[alj(xj)2(axjb)lj(xj)]0.解出

a2lj(xj),b12xjlj(xj). 由于

lj(x)

xxj1xxj1xx0xx1xxn

,

xjx0xjx1xjxj1xjxj1xjxn

n

故

lj(xj)

k0kj

1

, xjxk

n

于是

j(x)[12(xxj)

k0kj

1

]l2j(x).xjxk

(5)

同理，可得

j(x)(xxj)l2j(x). (6) 将(5)式、(6)式 代入式(4)便得到埃尔米特插值多项式

n

12

H2n1(x)yi[12(xxi)]lj(x)yj(xxj)l2j(x)

xxj0k0jj0k

kj

n

n

(7)

满足条件(1)的埃尔米特插值多项式是唯一的。这可用反证法证明，此处

从略。

2、两点三次埃尔米特插值

设已知yf(x)在[a,b]上的节点x0,x1上的函数值y0,y1及一阶导数值

,y1,则可按公式(7)写出三次埃尔米特插值多项式 y0

0(x)y11(x)H3(x)y00(x)y11(x)y0 y0(12

xx0xx12xx1xx02

)()y1(12)() x0x1x0x1x1x0x1x0

1

y(xx)(xx

x0x1

(xx1)()2y1

xx02

).x1x0

二、算法框图

两点三次埃尔米特插值框图如下：

两点三次埃尔米特插值框图如下：

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.[**************]