数值分析 第2章 插值法

Ch2插值法2.1引言2.2拉格朗日插值2.3均差与牛顿插值公式2.4差分与等距节点插值2.5埃尔米特插值2.6分段低次插值2.7三次样条插值当精确函数y=f(x)非常复杂或未知时，在一系列节点x0…xn

处测得函数值y0

=f(x0),…yn

=f(xn)，由此构造一个简单易算的近似函数p(x)

f(x)，满足条件p(xi)=f(xi)(i=0,…n)。这里的p(x)

称为f(x)的插值函数。最常用的插值函数是…?多项式§2.1引言x0x1x2x3x4xg(x)

f(x)点称为插值节点.求插值函数p(x)的方法称为插值法.1.若p(x)是次数不超过n的代数多项式,即为插值多项式.2.若p(x)为分段的多项式,就称为分段插值.3.若p(x)为三角多项式,就称为三角插值.√√插值多项式就是根据给定的n+1个点求一个n次多项式:使得：p(xi)=yi(i=0,1,…,n)即:这里是n+1个待定系数,根据n+1个条件得到的方程组是关于参数的线性方程组.当节点互异时由于系数行列式所以解是唯一存在的,但直接求解较复杂,也得不到统一的表达式.所以通常求插值多项式不用这种方法.而使用下节给出的基函数方法.返回主页§2.2拉格朗日插值

/*LagrangePolynomial*/P（xi)=yi,i=0,1,…,n条件：无重合节点，即n=1已知x0

,x1

;

y0

,

y1

，求使得111001)(,)(yxLyxL==可见L1(x)是过(x0,y0

)和(x1,y1

)两点的直线。)()(0010101xxxxyyyxL---+=101xxxx--010xxxx--=y0

+y1l0(x)l1(x)==10)(iiiyxl称为拉氏基函数

/*LagrangeBasis*/，满足条件li(xj)=ij

/*KroneckerDelta*/求n

次多项式使得l0(x)及l1(x)为线性插值基函数,满足l0(x0)=1,l0(x1)=0l1(x0)=0,l1(x1)=1n=2假定插值节点为,要求二次插值多项式L2(x),满足L2(xj)=yj

(j=0,1,2)可采用基函数方法,此时基函数l0(x),l1(x),l2(x)是二次函数,且在节点满足l0(x0)=1,l0(x1)=0,l0(x2)=0l1(x0)=0,l1(x1)=1,l1(x1)=0l2(x0)=0,l2(x1)=0,l2(x2)=1满足上述条件的基函数是很容易求出的.例如求l0(x),因为它有两个零点x1,x2,故可表示为l0(x)=A(x-x1)(x-x2）,由l0(x0)=1,可以解出于是:同理所以二次插值多项式为:n

2希望找到li(x)，i=0,…,n

使得

li(xj)=ij

；然后令==niiinyxlxL0)()(，则显然有Ln(xi)=

yi

。li(x)每个li有n

个根x0…

xi…xn=-=---=njjijiniiixxCxxxxxxCxl00)())...()...(()(-==jijiiiixxCxl)(11)(LagrangePolynomial引入记号则定理(唯一性)满足的n

阶插值多项式是唯一存在的。证明：反证：若不唯一，则除了Ln(x)外还有另一n

阶多项式Pn(x)满足Pn(xi)=yi

。考察则Qn

的阶数n而Qn有个不同的根n+1x0…xn注：若不将多项式次数限制为n

，则插值多项式不唯一。例如也是一个插值多项式，其中可以是任意多项式。

插值余项与误差估计定理2设在[a,b]上连续,在(a,b)内存在,节点是满足条件的插值多项式,则对任何插值余项这里且依赖于x,设节点在[a,b]内存在,考察截断误差，且f

满足条件,Rolle’sTheorem:若充分光滑，，则存在使得。推广：若使得使得Rn(x)至少有个根n+1=-=niinxxxKxR0)()()(任意固定x

xi(i=0,…,n),考察=-=niixtxKtRnt0)()()()(j(x)有n+2

个不同的根x0…

xn

x!)1()()()1(+-+nxKRxnnx注意这里是对t求导=+--++!)1)(()()()1()1(nxKLfxnnxnxx!)1()()()1(+=+nfxKxnx注：

通常不能确定x

,而是估计,x(a,b)

将作为误差估计上限。当

f(x)为任一个次数n

的多项式时，,可知，即插值多项式对于次数n的多项式是精确的。Quiz:

给定xi=i+1,i=0,1,2,3,4,5.

下面哪个是l2(x)的图像？

y

0

-

-

-

1

0.5

-0.5

1

2

3

4

5

6

x

y

0

-

-

-

1

0.5

-0.5

1

2

3

4

5

6

x

y

0

-

-

-

1

0.5

-0.5

1

2

3

4

5

6

x

ABC例：已知分别利用sinx的1次、2次Lagrange插值计算sin50

并估计误差。解：n=1分别利用x0,x1

以及x1,x2

计算利用这里而sin50=0.7660444…)185(50sin10pL0.77614实际误差0.01001利用sin500.76008,实际误差0.00596n=2)185(50sin20pL0.76543sin50=0.7660444…2次插值的实际误差0.00061注：1.高次插值通常优于低次插值；2.但绝对不是次数越高就越好。例题：已知函数表X1.12751.15031.1735F(X)0.11910.139540.15932应用1次、2次Lagrange计算F(1.1300)的近似值。解：二次插值functiony=lagrange(x0,y0,x)n=length(x0);m=length(x);fori=1:mz=x(i);s=0.0;fork=1:np=1.0;forj=1:nifj~=kp=p*(z-x0(j))/(x0(k)-x0(j));endends=p*y0(k)+s;endy(i)=s;end拉格朗日插值函数：x=[pi/6,pi/4,pi/3];>>y=[1/2,1/sqrt(2),sqrt(3)/2];>>x0=5*pi/18;>>y0=lagrange(x,y,x0)y0=0.76543389522903返回主页x1=[pi/6,pi/4];>>y1=[1/2,1/sqrt(2)];>>x0=5*pi/18;>>y11=lagrange(x1,y1,x0)y11=0.77614237491540>>x2=[pi/4,pi/3];>>y2=[1/sqrt(2),sqrt(3)/2];>>y12=lagrange(x2,y2,x0)y12=0.76007965538584上述三种情形在MATLAB的演示：§2.3均差与牛顿插值公式

/*Newton’sInterpolation*/Lagrange插值虽然易算，但若要增加一个节点时，全部基函数li(x)都需重新算过。的形式，希望每加一个节点时，将Ln(x)改写成只附加一项上去即可。????其中均为待定系数,可由插值条件确定.Pn(x)当x=x0时,当x=x1时,当x=x2时,依次可递推得到为写出系数ak的一般表达式,先引进均差的定义.均差(亦称差商)

/*divideddifference*/定义2

称为函数f(x)关于点x0,xk的一阶均差.为函数f(x)的二阶均差.一般地,称为函数f(x)的k阶均差./*the1stdivideddifferenceoffw.r.t.x0

andxk

*/均差性质1:k阶均差可表示为函数值f(x0),…,f(xk)的线性组合,即(3.3)注:这性质表明均差与节点的排列次序无关,称为均差的对称性.(3.2)当k=1时,假设当k=n时,结论成立.当k=n+1时,用归纳法进行证明证毕.均差性质2由性质1及(3.2)可得(3.4)证明:均差性质3若f(x)在[a,b]上存在n阶导数,且节点则n阶均差与导数关系如下:(3.5)证明:设处均为0,所以在[a,b]上有n+1个零点,根据罗尔定理,在[a,b]内至少有n个零点.反复应用罗尔定理,可知在[a,b]内至少有一个零点.记为使均差计算可列均差表,如下xkf(xk)一阶均差二阶均差三阶均差四阶均差x0f(x0)x1f(x1)f[x0，x1]x2f(x2)f[x1，x2]f[x0，x1，x2]x3f(x3)f[x2，x3]f[x1，x2，x3]f[x0，x1，x2,x3]x4f(x4)f[x3，x4]f[

x2，x3,x4]f[x1，x2，x3,x4]f[x0，x1，x2,x3，x4]………………牛顿插值

/*Newton’sInterpolation*/12…………n+11+(x

x0)2+……+(x

x0)…(x

xn1)n+1Nn(x)Rn(x)ai=

f[x0,…,xi]满足插值条件，且次序不超过n．注：

由唯一性可知Nn(x)Ln(x)，只是算法不同，故其余项也相同，即实际计算过程为:f(x0)f(x1)f(x2)…f(xn1)f(xn)f[x0,x1]f[x1,x2]…………f[xn1,xn]f[x0,x1,x2]…………f[xn2,xn1,xn]f[x0,…,xn]

f(xn+1)f[xn,xn+1]f[xn1,xn,xn+1]f[x1,…,xn+1]f[x0,…,xn+1]例２．给出f(x)的函数表,求4次牛顿插值多项式,并由此计算f(0.596)的近似值.解:首先根据给定的函数表构造出均差表:xkf(xk)一阶均差二阶均差三阶均差四阶均差0.400.410750.550.578151.116000.650.696751.186000.280000.800.888111.275730.358730.197330.901.026521.384100.433480.213000.031341.051.253821.515330.524930.228630.03126-0.00012N4(x)=0.40175+1.116(x-0.4)+0.28(x-0.4)(x-0.55)+0.19733(x-0.4)(x-0.55)+0.19733(x-0.4)(x-0.55)(x-0.65)+0.03134(x-0.4)(x-0.55)(x-0.65)(x-0.8)于是,f(0.596)≈N4(0.596)=0.63192截断误差|R4(x)|≈|f[x0,…,x5]w5(0.596)|≤3.63*10-9functionf=newton(x,y,a)m=length(x);n=length(a);fori=1:m-1q(i,1)=(y(i+1)-y(i))/(x(i+1)-x(i));endforj=2:m-1fori=j:m-1q(i,j)=(q(i,j-1)-q(i-1,j-1))/(x(i+1)-x(i+1-j));endendb(1)=y(1);fori=2:mb(i)=q(i-1,i-1);endfori=1:ns=0;forj=1:mk=b(j);forl=1:j-1k=k*(a(i)-x(l));ends=s+k;endf(i)=s;end>>x=[0.40.550.650.80.91.05];>>y=[0.410750.578150.696750.888111.026521.25382];>>x0=0.596;>>y0=newton(x,y,x0)y0=0.63191749923175返回主页在MATLAB演示：§2.4差分与等距节点公式向前差分/*forwarddifference*/kffk+1fk-=kkkkkkkkffff1111)(-+---==向后差分/*backwarddifference*/111----=kkkkkkfffk1kkfff-=中心差分/*centereddifference*/其中设函数y=f(x)在等距节点xk=x0+kh(k=0,1,…,n)上的值fk=f(xk)为已知,这里h为常数,成为步长.向前差分算子；向后差分算子；δ中心差分算子此外还可定义不变算子I及位移算子E为：同理可得于是，由

,可得

由差分定义并应用算子符号运算可得下列基本性质.性质1各阶差分均可用函数值表示.例如

其中

为二项式展开系数性质2

可用各阶差分表示函数值.例如，可用向前差分表示因为于是性质3

均差与差分有的关系.由定义可知，向前差分一般地有同理，对向后差分有(4.7)(4.8)利用(4.7)及(3.5)又可得到(4.9)这就是差分与导数的关系.计算差分可列差分表fkf0f1f2f3f4……………

2.4.2等距节点插值公式牛顿均差插值多项式(3.6)中各阶均差用相应差分代替，就可得到各种形式的等距节点插值公式.这里只推导常用的向前插值与向后插值公式。如果有节点于是将此式及(4.7)代入(3.6)，则得称为Newton前插公式，其余项由(2.14)得(4.10)要计算附近点x的函数f(x)的值，

可令如果要用函数表示xn附近的函数值f(x)，此时应用牛顿插值公式(3.6)，插值点应按xn,xn-1,…,x0的次序排列，有

称为牛顿向后插值公式,其余项作变换并利用公式(4.8)，代入上式得

，

试用三次等距节点插值公式求f(1.01)及f(1.28)的近似值.

例3

设给出在的值.解本题只要构造出f的差分表，再按Newton向前插值公式及向后插值公式计算即可.的差分表如图所示.计算f(1.01)可用Newton向前插值公式(4.10)，此时用到差分表中的上半部分划波纹线的各阶差分值计算f(1.28)要用Newton向后插值公式(4.12)，它用到差分表下部分的差分Δ(即下划直线的).f(1.01)与f(1.28)的7位有效数字分别为可见计算结果已相当精确牛顿公式牛顿前差公式/*Newton’sforward-differenceformula*/牛顿后差公式/*Newton’sbackward-differenceformula*/将节点顺序倒置：设，则)()()(000xfkthtxNxNknknn=+==设，则)()1()()(0nknkknnnxfkthtxNxN--=+==注：一般当x

靠近x0时用前插，靠近xn时用后插，故两种公式亦称为表初公式和表末公式。返回主页§2.5hermite插值引入

在实际问题中，对所构造的插值多项式，不仅要求函数值重合，而且要求若干阶导数也重合。即：要求插值函数P(x)满足:（1）把此类插值多项式称为埃米尔特（Hermite）插值多项式或称带导数的插值多项式，记为H(x)。H(x)

存在且唯一推导只讨论函数值与导数值个数相等的情况.设在节点上,要求插值多项式,满足条件(2)这里给出的个条件,可唯一确定一个次数不超过的多项式其形式为如根据条件(2)来确定个系数,显然非常复杂,因此,仍采用求Lagrange插值多项式的基函数方法.先求插值基函数及,共有个,每一个基函数都是次多项式,且满足条件(3)于是满足条件(2)的插值多项式可写成用插值基函数表示的形式,即(4)由条件(3),显然有下面的问题就是求满足条件(36)的基函数及为此,可利用Lagrange插值基函数.令其中是式

所表示的基函数.由条件式(3),有整理,得解得由于两端取对数再求导,得于是同理可得(5)(6)

Hermite插值余项还可以证明满足条件(2)的插值多项式是唯一的.用反证法,假设及均满足条件(2),于是在每个节点上均有二重根,即有重根.但是不高于次的多项式,故

.唯一性得证.(7)其中且与有关.仿照Lagrange插值余项的证明方法,若f(x)在(a,b)内的2n+1阶导数存在，则其插值余项为多少？作为带导数插值多项式(4)的重要特例是n=1的情形.这时可取节点为及,插值多项式为,满足条件(8)相应的插值基函数为,它们满足条件根据(5)式及(6)式的一般表达式,可得于是满足条件(8)的插值多项式是其余项为,由式(7)得注：

N

个条件可以确定阶多项式。其余项为N

1要求在1个节点x0处直到m0阶导数都重合的插值多项式即为在点处的Taylor多项式举例例：设

x0

x1x2,已知

,

求多项式

P(x)满足

解：首先，P的阶数为3,模仿Lagrange多项式的思想，设其中且

,并估计误差。h0(x)有根x1,x2，且

x1是重根。又:h0(x0)=1C0

h2(x)与h0(x)完全类似。h1(x)有根

x0,x2

由余下条件

h1(x1)=1和

可解。

(x)h1

有根

x0,x1,x2又:C1可解。

与Lagrange分析完全类似解:设hi(x)有根x0

,…,xi,…,xn且都是2重根一般地，已知x0

,…,xn

处有y0

,…,yn

和，求

H2n+1(x),满足

。其中由余下条件和

可解Ai

和Bi

(x)hi有根x0

,…,xn,除了xi

外都是2重根设则EXE

给定试求f(x)在[1/4,9/4]上的三次Hermit插值多项式P(x),使它满足并写出余项表达式.Quiz:

给定xi=i+1,i=0,1,2,3,4,5.

下面哪个是h2(x)的图像？x0--10.5123456yxy0---10.5123456斜率=1

求Hermite多项式的基本步骤：写出相应于条件的hi(x)、hi(x)的组合形式；对每一个hi(x)、hi(x)找出尽可能多的条件给出的根；根据多项式的总阶数和根的个数写出表达式；根据尚未利用的条件解出表达式中的待定系数；最后完整写出H(x)。返回主页§2.6分段低次插值

/*piecewisepolynomialapproximation*/多项式插值的龙格现象例：在[5,5]上考察的Ln(x)。取

-5

-4

-3

-2

-1

0

1

2

3

4

5

-0.5

0

0.5

1

1.5

2

2.5

Ln(x)f(x)n

越大，端点附近抖动越大，称为Runge现象返回主页分段线性插值

/*piecewiselinearinterpolation*/在每个区间上，用1阶多项式

(直线)逼近f(x):记，易证：当时，一致失去了原函数的光滑性。分段Hermite插值

/*Hermitepiecewisepolynomials*/给定在上利用两点的y及y’构造3次Hermite函数导数一般不易得到。返回主页§2.7三次样条插值

/*CubicSpline*/

样条(spline),是早期飞机、造船工业中绘图员用来画光滑曲线的细木条或细金属丝.绘图时,为了将一些已知的点连成光滑曲线,绘图员用压铁把样条固定在相邻的若干点上,样条具有弹性,形成通过这些点的光滑曲线,沿着它就可画出所需的曲线.数学上仿此得出的函数便称为样条函数.

样条插值是用分段低次多项式去逼近被逼近函数,并且能够满足对光滑性的要求,又无需给出每个节点处的导数值.它除了要求给出各点处的函数值之外,只需提供两个边界节点处的导数信息.引入要求：插值曲线既要简单，又要在曲线的连接处比较光滑。

这样的分段插值函数在分段上要求多项式次数低，这种插值方法称为——样条插值。它所对应的曲线称为样条曲线，其节点称为样点，把满足这样条件的插值函数，称为样条插值函数，而在节点上不仅连续，还存在连续的低阶导数，定义设。三次样条函数,

且在每个上为三次多项式

/*cubicpolynomial*/。若它同时还满足，则称为f的三次样条插值函数

/*cubicsplineinterpolant*/.注：三次样条与分段Hermite插值的根本区别在于S(x)自身光滑，不需要知道f的导数值（除了在2个端点可能需要）；而Hermite插值依赖于f在所有插值点的导数值。f(x)H(x)S(x)

从定义知,S(x)在每个小区间上要确定4个待定系数,共有n个小区间,故应确定4n个参数.根据S(x)在[a,b]上二阶导数连续,在节点xj(j=1,2,…,n-1)处应满足连续性条件

S(xj-0)=S(xj+0),S’(xj-0)=S’(xj+0),S’’(xj-0)=S’’(xj+0),

共3n-3个条件,再加上S(x)满足插值条件S(xj)=yj(j=0,1,…,n),共n+1个,共有4n-2个条件,

因此还需要2个条件才能确定S(x).通常可在区间[a,b]的端点各加一个条件(称为边界条件).分析：构造三次样条插值函数的三弯矩法

/*methodofbendingmoment*/在上，记],[for)()(1][jjjxxxxSxS-=对每个j,此为3次多项式则S[j]”(x)为次多项式，需个点的值确定之。12设S[j]”(xj1)=Mj1，S[j]”(xj)=Mj

对应力学中的梁弯矩，故名对于x

[xj1,

xj

]可得到S[j]”(x)=jjjjjjhxxMhxxM11---+-积分2次，可得S[j]’(x)和S[j](x)：jjjjjjjAhxxMhxxM+-+----2)(2)(2121S[j]’(x)=jjjjjjjjBxAhxxMhxxM++-+---6)(6)(3131S[j](x)=利用已知S[j](xj1)=yj1S[j](xj)=yj

可解jjjjjjjhMMhyyA611-----=jjjjjjjjjjjjhxxhMyhxxhMyBxA12211)6()6(-----+--=+jjjjjjhxxMhxxM+-+---6)(6)(3131S[j](x)=jjjjjjjjjjhxxhMyhxxhMy12211)6()6(-----+--类似地可以求出S(x)在区间

在上，记],[for)()(j+1][jxj+1xxxSxS=则S[j+1]”(x)为次多项式，需个点的值确定之。12设S[j+1]”(xj)=Mj，S[j+1]”(xj+1)=Mj+1

对于x

[xj,

xj+1

]可得到S[j+1]”(x)=J+1jJ+1J+1J+1jhxxMhxxM-+-积分2次，可得S[j+1]’(x)和S[j+1](x)：jJ+1jJ+1J+1jAhxxMj+1hxxM+-+--2)(2)(22S[j+1]’(x)=jjJ+1jJ+1J+1J+1jBxAhxxMhxxM++-+-6)(6)(33S[j+1](x)=利用已知S[j+1](xj)=yj

S[j+1](xj+1)=yj+1

可解的表达式下面解决Mj

：利用S’

在xj的连续性[xj1,

xj

]:

S[j]’(x)=jjjjjjjjjjjhMMxxfhxxMhxxM6],[2)(2)(112121------+-+--1111211216],[2)(2)(+++++++--+-+--jjjjjjjjjjjhMMxxfhxxMhxxM[xj,

xj+1]:

S[j+1]’(x)=利用S[j]’(xj-0)=S[j+1]’(xj+0)，合并关于Mj1、

Mj、Mj+1的同类项，并记,,,整理后得到：11jjjjhhh+++=l1jj-=lm]),[],[(6111jjjjjjjxxfxxfhhg-++-+=211gMMMjjjjjj=+++-lm

j

1n1即：有个未知数，个方程。n1n+1还需2个边界条件

/*boundaryconditions*/

第1类边条件

/*clampedboundary*/：S’(a)=y0’，S’(b)=yn’[a,

x1

]:

S[1]’(x)=1011012112106],[2)(2)(hMMxxfhaxMhxxM--+-+--010110)],[(62gy0’xxfhMM=-=+nnnnnngxxfyn’hMM=-=+--]),[(6211类似地利用[xn1,

b]上的

S[n]’(x)令：矩阵表示为:

第2类边条件：S”(a)=y0”=

M0，S”(b)=yn”=

Mn这时：特别地，M0=

Mn=0

称为自由边界

/*freeboundary*/，对应的样条函数称为自然样条

/*NaturalSpline*/。

第3类边条件/*periodicboundary*/

：当f

为周期函数时，

yn

=y0

，

S’(a+)=S’(b)

M0=

Mn

误差与收敛性定理4

构造三次样条插值函数的三转角算法下面构造用一阶导数值表示的三次样条插值函数.在力学上解释为细梁在截面处的转角,并且得到的转角与相邻两个转角有关,故称该算法为三转角算法.二次导数得:于是有:可选一种边界条件补充2个方程.对于边界条件1,则对于边界条件2,则可导出两个方程对于边界条件3,可得>>a=[01491625364964];>>b=[012345678];>>x=0:0.1:64;>>y=spline3(a,b,x,1,0,1/16);>>plot(x,y,'r');>>holdon;>>plot(a,b,'b');>>xla