牛顿插值和拉格朗日插值

数值计算方法

Newton插值函数和三次插值样条函数

主要内容

一、插值函数的基本思想

二、插值函数的基本算法

三、插值函数的Matlab实现

四、举例说明

五、插值函数的优、缺点

六、附录

一、插值函数的基本思想

在科学研究和实践活动中，我们所遇到的大量函数，有相当一部分是通过观测和实验得到的，虽然其函数关系在某个区间上是存在的，但是却不知道具体的表达式，只能通过观测或实验得到一些离散点的函数值、导数值，因此希望对这样的函数用一个简单的解析表达式近似的给出整体上的描述。

在数学上，该问题相当于，函数()x f y =在某个区间[]b a ,上存在且连续，但不知道解析表达式，只给出[]b a ,上离散点i x 的函数值

()i x f y =，如何计算()x f y =在[]b a ,上其它点的函数值。

插值法的基本思想就是构造一个简单函数()x P y =作为()x f 的近似表达式，以()x P 的值作为函数()x f 的近似值，而且要求()x P 在给定点i x 与取值相同，即()()i i x f x P =。通常我们称()x P 为()x f 的插值函数，其中i x 称为插值节。

二、插值函数的基本算法

1、Newton 插值函数 1）Newton 插值公式

首先考察待定常数()n i c i ,,2,1 =与差商之间的关系。

[]101,x x f c =，[]2102,,x x x f c =。

对于一般的()n j c j ,,4,3 =，我们用归纳法证得

[]

j j x x f c ,,0 =。

因此Newton 插值公式 为

()()[]()[]()()()

110100100,,,----++-+=n n n x x x x x x x x x f x x x x f x f x N

。

2）Newton 向前插值公式

设插值节点为等距节点n i ih x x i ，，

，100=+=，在newton 插值公式中用向前差分代替差商，可得

()()()()()()()()1100

1020

2000!!2!1-----?+--?+-?+=n i n

n n x x x x x x x x h

n f x x x x h

f x x h f f x N

令[]n t th x x ，，00∈+=，则()n i h i t x x i ,,1,0 =-=-，，于是

()()

()()().

11!

1!2!10

02000---?+-?+?+=+=n t t t n f t t f t f f th x N x N n n n 3）Newton 向后插值公式

在Newton 插值公式中改变插值节点的次序并用向后差分代替差商，可得

()[]()[]()()()

()()()()()()1112211011!!2!1,,,x x x x x x h

n f x x x x h f x x h f f x x x x x x x x x f x x x x f f x N n n n

n

n n n n

n n n n n n n n n n n n ---?++--?+-?+=---+-+=-----

令th x x n +=，则可得

()()

()()().

11!

1!2!12-++?+++?+?+=+=n t t t n f t t f t f f th x N x N n

n n n n n n n 注：① 差商的定义

称

[][][]

i

k i k i i k i i k i i i x x x x f x x f x x x f --=

+-+++++111,,,,,

为()x f 关于点k i i i x x x ++ ,,1的k 阶差商，其中称[]()j j x f x f =为()x f 关于j x 的零阶差商。 ② 差分的定义

设函数()x f y =在等距节点()n k kh x x k ,2,10=+=上的函数值

()i i x f f =，()

n

a b h -=

称为步长。

() 2102

212121±±=??? ??==±=±±±，，

，，，i x f f x f f h

x x i i i i i i 引入记号

2

1211

1-+-+-=-=?-=?i i i i i i i i i f f f f f f f f f δ

分别称()i x f 为在节点i x 的一阶向前差分，一阶向后差分和一阶中心差分，δ，，??称为差分算子。

2、三次插值样条函数（三弯矩方程） 首先引入区间[]b a ,的一个划分

b x x x x a n n =<<<<=?-110: ，

称n n x x x x ,,,110- 为剖分点，11,-n x x 为内剖分点，小区间[]j j j x x e ,1-=为

第j 个剖分单元。

记()()()i i i i i i x f f x f f x f f ''='''='=,,，它们分别表示被插函数()x f 在节点i x 处的函数值、一阶和二阶导数值。

设关于剖分?的分段多项式插值函数空间为

(){

[]}n i e x P s b a C s l k Sp i k ,2,1,,:,;;1=∈∈∈=?，

它表示在每个剖分单元为k 次多项式且在整个定义区间[]b a ,上l 阶连续可微的函数的全体。称()?;;l k Sp 为关于剖分?重数为l k -的k 次多项式样条空间。其空间维数为()()()111-+-+n l k n 。

这时取插值函数空间为()?;3Sp ，它表示在每个剖分单元为三次多项式且在区间[]b a ,上二阶连续可微的函数全体。易知该插值函数空间的维数为()3134+=--n n n 。

利用该空间，可以定义相应的三次插值样条函数()?∈;3Sp s ，其关键是给出适当的插值条件，首先自然要求()x s 在所有内节点上满足

()1,2,1,-==n i f x s i i

但这些插值条件的个数仅为1-n ，为了将()x s 唯一确定，还需要提供另外4个插值条件（这里我们只考虑第一种取法） 即I 型边界条件：

()()()()'

=''='==n n n n f x s f x s f x s f x s ,;,0000

通过计算，若记

1

1+++=

j j j j h h h λ，1

1++=

-=j j j j j h h h λμ，1,2,1-=n j

我们有

j j j j j j d M M M =+++-112λμ，1,2,1-=n j

这里[]11,,6+-=j j j j x x x f d ，其中[]

11,,+-j j j x x x f 表示()x f 的二阶差商。

在端点0x 和n x 处，我们利用附加边界条件()()n n f x s f x s '=''='1001

,，有

???

????=???? ??--'=+=???

? ??'--=+--n n n n n n n n d h f f f h M M d f h f f h M M :62:62110010

1110

从而可以得到一个关于待定参数n M M M ,,10的1+n 阶线性代数方程组，其矩阵形式为

???

???

???

???????????=??????????????

?????????????????????????

?----n n n n n n d d d d d M M M M M 121011011

2

21

122122

2

12

λμλμλμ（*） 我们称方程组（*）为三弯矩方程。

三、插值函数的

Matlab 实现

1）Newton 插值函数

function []=newton(x,y,v)

x=input('X='); y=input('Y='); v=input('V='); n=length(x); t=zeros(n,n); u=0;

for i=1:n

t(i,1)=y(i); end

for j=2:n for i=2:n if i>=j

t(i,j)=(t(i,j-1)-t(i-1,j-1))/(x(i)-x(i-j+1)); end end end

for k=1:n s=1;

for j=1:k

if j

s=s*(v-x(j));

end

end

m=s*t(k,k);

u=u+m;

end

disp('插值结果=');

disp(u);

end

2）Newton向前插值函数

function y0=newtonint1(x,y,x0)

n=length(y);

h=x(2)-x(1);

t=(x0-x(1))/h;

Y=zeros(n);Y(:,1)=y';

for k=1:n-1

Y(:,k+1)=[diff(y',k);zeros(k,1)]; end

y0=Y(1,1);

for i=1:n-1

z=t;

for k=1:i-1

z=z*(t-k);

end

y0=y0+Y(1,i+1)*z/prod([1:i]);

end

3）Newton向后插值函数

function y0=newtonint2(x,y,x0)

n=length(y);

h=x(n)-x(n-1);

t=(x(n)-x0)/h;

Y=zeros(n);

Y(:,1)=y';

for k=1:n-1

Y(:,k+1)=[zeros(k,1);diff(y',k)]; end

y0=Y(n,1);

for i=1:n-1

for k=1:i-1

z=z*(t-k);

end

y0=y0+Y(n,i+1)*(-1)^i*z/prod([1:i]);

end

4）三次插值样条函数（三弯矩方程）

function y0=csyt(x,y,z,x0)

n=length(x);m=length(y);

if n~=m

error;

end

if isempty(z)==1

z=[(y(2)-y(1))]/(x(2)-x(1))

(y(n)-y(n-1))/(x(n)-x(n-1));

end

h=zeros(1,n);l=ones(1,n);m=ones(1,n);M=zeros(n,1);d=zeros (n,1);

for k=2:n

h(k)=x(k)-x(k-1);

end

for k=2:n-1

l(k)=h(k+1)/(h(k)+h(k+1));m(k)=1-l(k);

d(k)=6/(h(k)+h(k+1))*((y(k+1)-y(k))/h(k+1)-(y(k)-y(k-1))/ h(k));

end

d(1)=6/h(2)*((y(2)-y(1))/h(2)-z(1));

d(n)=6/h(n)*(z(2)-(y(n)-y(n-1))/h(n));A=diag(2*ones(1,n)) ;

for i=1:n-2

A(i,i+1)=l(i);A(i+1,i)=m(i+1);

M=A\d;

end

for k=2:n

if x(k-1)<=x0*x0<=x(k)

y0=M(k-1)/6/h(k)*(x(k)-x0)^3+M(k)/6/h(k)...

*(x0-x(k-1))^3+1/h(k)*(y(k)-M(k)*h(k)^2/6)...

*(x0-x(k-1))+1/h(k)*(y(k-1)-M(k-1)*h(k)^2/6)*(x(k)-x0);

return

end

end

四、举例说明

例1 已知函数()()x sh x f =的函数值如下表所示，

i x

0.4 0.55 0.65 0.80 0.90 ()i x f

0.41075

0.57815

0.69675

0.88811

1.02652

构造4次Newton 插值多项式并计算()596.0f 的值。

解 在Matlab 命令窗口执行 >> newton

X=[0.40 0.55 0.65 0.80 0.90];

Y=[0.41075 0.57815 0.69675 0.88811 1.02652]; V=0.596 插值结果= 0.6319

另外，如果我们动手自己算，则4次Newton 插值多项式为

()()()()()()()()()()()80.065.055.04.003134

.065.055.04.019733

.055.04.028000.04.011600

.141075.04----+---+--+-+=x x x x x x x x x x x N

于是

()()63195.0596.0596.04=≈N f 。

其中构造差商表如下：

i x

()i x f

一阶差商 二阶插商 三阶差商 四阶差商

0.40 0.41075 0.55 0.57815 1.11600 0.65 0.69675 1.18600 0.28000 0.80 0.88811 1.27573 0.35893 0.19733 0.90 1.02652

1.38410

0.43348

0.21300

0.03134

例2 已知函数()x f 的函数表如下:

i x

1 2 3 4 5 6

i y

1.0000 1.2599 1.4422 1.5874 1.7100 1.8171

求函数()x f 在6.5=x 处的函数值。

解 ①用Newton 向前插值公式，在Matlab 命令窗口执行 >> x=[1 2 3 4 5 6];

>> y=[1.0000 1.2599 1.4422 1.5874 1.7100 1.8171]; >> x0=5.6;

>> y0=newtonint1(x,y,x0) y0 =

1.7754

解 ②用Newton 向后插值公式，在Matlab 命令窗口执行 >> x=[1 2 3 4 5 6];

>> y=[1.0000 1.2599 1.4422 1.5874 1.7100 1.8171]; >> x0=5.6;

>> y0=newtonint2(x,y,x0)

y0 =

1.7754

例3 已知()x f 在16941，，，

=x 点的函数值为()4321，，，=x f ，试求函数()x f 满足边界条件()()58

1

16211=='='x s s ，，的三次样条。

解 在Matlab 命令窗口执行 >> x=[1 4 9 16]; >> y=[1 2 3 4]; >> z=[1/2 1/8]; >> x0=5;

>> y0=csyt(x,y,z,x0) y0 =

2.2366

四、插值函数的优缺点 1、Newton 插值函数

优点：从Newton 插值公式可以看出，当增加一个节点时，公式只是在原来的基础上增加一项，因此前面计算的结果仍然可以用到。即具有“承袭性”。所以Newton 插值函数具有灵活增加节点的优点。

缺点：根据直观想象和插值余项估计，似乎插值节点的个数越多，插值多项式次数越高，与被插函数越接近。事实并非如此。实际应用时，高次插值（如7、8次以上）的逼近效果并不好，节点的增多当然能使被插函数()x P 在更多地方与()x f 相等，但是在两个节点之间

()x P 不一定能很好的逼近()x f ，有时差别很大。例如龙格现象（备注）。

2、三次插值样条函数

优点：高次插值函数的计算量大,有剧烈振荡,且数值稳定性差；在分段插值中，分段线性插值在分段点上仅连续而不可导，分段三次Hermite 特插值有连续的一阶导数，如此光滑程度常不能满足物理问题的需要，三次插值样条函数可以同时解决这两个问题，使插值函数既是低阶分段函数，又是光滑的函数，并且只需在区间端点提供某些导数信息。

缺点：三次插值样条函数公式不能局部确定，它需要求解一个线性代数方程组。

六、附录

备注：龙格函数

()2

11

x x f +=

，[]55，

-∈x Matlab 程序如下：

function a=Runge(n,x)

X=linspace(-5,5,n+1)';Y=1./(1+(X.*X));a=0; for i=1:n+1 li=1;

for j=1:n+1 if i~=j

li=li*(x-X(j))/(X(i)-X(j)); end end

a=a+li*Y(i); End

在Matlab 命令窗口执行：

>>X=linspace(-5,5,100)';Y=1./(1+X.*X); >> plot(X,Y,'r','LineWidth',1); >> hold on;

>> fplot('Runge(10,x)',[-5 5],'g'); >> fplot('Runge(5,x)',[-5 5],'b'); >> hold off; >> grid on