数值分析 Ch2-插值法

第二章插值法/*Interpolation*/Interpolation_introduction§1引言Interpolation_introduction§1引言1.函数表达式过于复杂不便于计算,而又需要计算许多点处的函数值2.仅有几个采样点处的函数值,而又需要知道非采样点处的函数值……上述问题的一种解决思路：建立复杂函数或者未知函数的一个便于计算的近似表达式.解决方法－插值法

Interpolation_introduction函数逼近的方法有很多，例如Taylor级数，Fourier级数，有限元方法、边界元方法，小波分析等，大学科叫逼近论。本书讨论连续函数的逼近，主要介绍插值法(chapter2)和最佳一直逼近、最小平方逼近离散数据拟合(chapter3)1.插值概念求插值函数(x)的问题(方法)称为插值问题(方法)。2.几何意义、内插法、外插法内插外插§2拉格朗日多项式

/*LagrangePolynomial*/niyxPiin,...,0,)(==求n次多项式使得条件：无重合节点，即l0(x)l1(x)§2LagrangePolynomial多项式插值是数值分析的基本工具，常用来计算被插函数的近似函数值，零、极点，导数、积分（第四章数值积分和数值微分），解微分方程（第五章）、积分方程x0x1x2x3x4xPn(x)

f(x)解几何上看，即求多项式曲线与被插值函数曲线间满足：特点：插值曲线Pn(x)过被插值曲线f(x)的上给定的n+1个点。n=1已知x0,x1;

y0,

y1，求使得111001)(,)(yxPyxP==可见P1(x)是过(x0,y0)和(x1,y1)两点的直线。)()(0010101xxxxyyyxP---+=101xxxx--010xxxx--=y0+y1l0(x)l1(x)==10)(iiiyxl称为拉氏基函数满足条件li(xj)=ij

§2LagrangePolynomial

提问：上面所提的多项式Pn(x)是否存在？若存在，是否唯一？如何求？证明插值条件(2.2)等价于线性方程组定理1

满足插值条件(2.2)的不超过n次的插值多项式唯一存在。系数行列式（n+1阶范德蒙行列式）由克莱默法则知，方程组有唯一解2-1插值多项式的存在唯一性

唯一性的另一证明满足的n阶插值多项式是唯一存在的。证明

(前面已利用Vandermonde行列式论证)反证：若不唯一，则除了Ln(x)外还有另一n阶多项式Pn(x)满足Pn(xi)=yi。考察则Qn的阶数n而Qn有个不同的根n+1x0…xn注：若不将多项式次数限制为n，则插值多项式不唯一。例如也是一个插值多项式，其中可以是任意多项式。§2LagrangePolynomialInterpolationpolynomial

2-2线性插值与抛物插值

x0x1(x0,y0)(x1,y1)P1(x)f(x)可见P1(x)是过(x0,y0

)和(x1,y1

)两点的直线。直线方程为:等价变形为:记为:1.线性插值引入记号:则:分析两个基函数有:对于三个点,类似有:称为插值基函数x0x1x2P2(x)

f(x)f(x)2.抛物线(二次)插值将以上思路推广到n+1个节点情形，即可得到类似的插值基函数和插值多项式表示形式。P2(x)xn

1希望找到li(x)，i=0,…,n使得li(xj)=ij；然后令，则显然有Pn(xi)=yi。li(x)每个li有n个根x0…

xi…xn==jiC0=-njijxx)(---inxxixxxxC0))...()...((ixl)(-=jijxixiC)(1=iixl1)(LagrangePolynomial与节点有关，而与f无关==niinxlxP0)()(yi基函数法（n=1情形的推广）§2LagrangePolynomial2-3Lagrange插值多项式

2-4插值余项

/*Remainder*/设节点在[a,b]内存在,考察截断误差，且f

满足条件,Rolle’sTheorem:若充分光滑，，则存在使得。推广：若使得使得存在使得Rn(x)至少有个根n+1=-=niinxxxKxR0)()()(任意固定xxi(i=0,…,n),考察=-=niixtxKtRnt0)()()()(j(x)有n+2个不同的根x0…

xn

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

通常不能确定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§1LagrangePolynomial例：已知分别利用sinx的1次、2次Lagrange插值计算sin50并估计误差。解：n=1分别利用x0,x1以及x1,x2计算利用这里而sin50=0.7660444…)185(50sin10pL0.77614外推

/*extrapolation*/的实际误差0.01001利用sin500.76008,内插

/*interpolation*/的实际误差0.00596内插通常优于外推。选择要计算的x所在的区间的端点，插值效果较好。§1LagrangePolynomialn=2)185(50sin20pL0.76543sin50=0.7660444…2次插值的实际误差0.00061高次插值通常优于低次插值但绝对不是次数越高就越好，嘿嘿……Whenyoustartwritingtheprogram,youwillfindhoweasyitistocalculatetheLagrangepolynomial.Ohyeah?WhatifIfindthecurrentinterpolationnotaccurateenough?Thenyoumightwanttotakemoreinterpolatingpointsintoaccount.Right.ThenalltheLagrangebasis,li(x),willhavetobere-calculated.Excellentpoint!Wewillcometodiscussthisproblemnexttime.§1LagrangePolynomial练习：设f(x)=x4,利用Lagrange插值写出以-1,0,1,2为节点的插值多项式。练习：已知插值节点满足xi-xi-1=h,i=1,2,3,证明三次插值多项式L3(x)与被插函数f(x)的差有如下关系：§1LagrangePolynomial练习：证明Interpolation_introduction问题的提出：如果需要增加精度，一般采用增加插值节点来实现。但是由于插值基函数的性质，使得计算新的插值多项式时，原来的计算量不能很好地利用，造成计算的浪费，为了克服这一缺点，我们将要介绍下面的逐步线性插值和Newton插值。§3逐次线性插值

/*LagrangePolynomial*/§3逐次线性插值法/*LagrangePolynomial*/实际上，是对两个低次插值的线性插值，这种通过低次插值再作线性插值生成高次插值的方法称为逐次线性插值。

Aitken法(按下表计算)线性插值基函数增加如果精度不够，增加节点x4,同时表中增加一行，三角形斜边上即为所要求的各次插值多项式。k1k0k2k3k4

Neville法(按下表计算)增加如果精度不够，增加节点x4,同时表中增加一行，三角形斜边上即为所要求的各次插值多项式。k1k0k1k1k1HW:用类似于前面的方法构造Neville计算公式注：Atkin方法和Neville方法与Lagrange公式相比，当需要增加节点时，很容易由低次插值构造高次插值，而Lagrange插值公式中，每个基函数都需要作适当变化。误差估计：由插值多项式的存在唯一性知，仍有但这里可采用一种更简便的方法。当f(n+1)(x)在插值区间变化不大时，设f(n+1)(x)L，则有根据前面的计算结果估计当前的误差：事后误差估计（实用），前面给出的误差估计（事先误差估计）不实用§4牛顿插值/*Newton’sInterpolation*/Lagrange插值虽然易算，但若要增加一个节点时，全部基函数li(x)都需重新算过。公式不具有继承性，不利于编程。将Ln(x)改写成的形式，希望加一个节点时，只附加一项上去即可。????

差商(亦称均差)

/*divideddifference*/1阶差商

/*the1stdivideddifferenceoffw.r.t.xi

andxj

*/2阶差商f(x0)1阶差商的几何意义：弦截线的斜率§4Newton’sInterpolation§4Newton’sInterpolation11101010111010],,...,[],,...,[],,...,[],...,,[],...,[++--+++--=--=kkkkkkkkkkkxxxxxfxxxfxxxxxfxxxfxxf(k+1)阶差商：事实上其中Warning:myheadisexploding…Whatisthepointofthisformula?差商的值与xi的顺序无关！1.线性：2.差商可以表示为函数值的线性组合：3.

对称性：由2知，差商的值与节点的顺序无关！4.差商的另一种定义：由2，3及均差定义可得§4Newton’sInterpolation§4Newton’sInterpolation牛顿插值

/*Newton’sInterpolation*/…………Nn(x)—n次多项式，满足：Nn(xi)=f(xi)Rn(x)—插值余项，满足Rn(xi)=0，i=0,…,n

ai=

f[x0,…,xi](1)(2)(n)(n)(n-1)…(2)(1)§4Newton’sInterpolation注：

由唯一性可知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]增加如果精度不够，增加节点xn+1,同时表中增加一行，三角形斜边上即为所要求的各次项系数。§1LagrangePolynomial练习：设f(x)=3x2+5,求f[x0,x1,x2],f[x0,x1,x2,x3]练习：练习：§4Newton’sInterpolation

等距节点公式

/*FormulaewithEqualSpacing*/向前差分

/*forwarddifference*/iiifff-=+1ikikikikffff1111)(-+---==向后差分

/*backwarddifference*/111----=ikikikfffi1iifff-=中心差分

/*centereddifference*/其中当节点等距分布时:fi=f(xi)

差分计算可通过构造差分表得到增加

差分的重要性质：

线性：例如

各阶差分可用函数值表示：其中/*binomialcoefficients*/§4Newton’sInterpolation

函数值可用各阶差分表示：

差商与差分的关系：

若f(x)是n次多项式，则是次多项式，而

差分与导数的关系(由差分与差商、差商与导数的关系得）：§4Newton’sInterpolation等距节点牛顿公式

牛顿前差公式/*Newton’sforward-differenceformula*/注：一般当x靠近x0时用前插，靠近xn时用后插，故两种公式亦称为表初公式和表末公式。

牛顿后差公式/*Newton’sforward-differenceformula*/例：§6埃尔米特插值/*HermiteInterpolation*/不仅要求函数值重合，而且要求若干阶导数也重合。即：要求插值函数(x)满足(xi)=f(xi),'(xi)=f'

(xi),…,(mi)(xi)=f

(mi)(xi).注：

n+1个条件可以确定指数不超过n次的多项式。要求在1个节点x0处直到m0

阶导数都重合的插值多项式即为Taylor多项式其余项为一般只考虑f与f'的值。埃尔米特插值§3HermiteInterpolation例：设x0

x1x2,已知f(x0)、f(x1)、f(x2)和f’(x1),求多项式P(x)满足P(xi)=f(xi)，i=0,1,2，且P’(x1)=f’(x1),并估计误差。模仿Lagrange多项式的思想，设解：首先，P

的阶数=3+=213)()()()()(=0iiixhx1f’xhxfxPh0(x)有根x1,x2，且h0’(x1)=0x1是重根。)()()(22100xxxxCxh--=又:h0(x0)=1C0h2(x)h1(x)有根x0,x2

))()(()(201xxxxBAxxh--+=由余下条件h1(x1)=1和

h1’(x1)=0可解。与h0(x)完全类似。

(x)h1有根x0,x1,x2

h1))()(()(2101xxxxxxCx---=h1又:’(x1)=1C1

可解。其中hi(xj)=ij,hi’(x1)=0,

(xi)=0,

’(x1)=1h1h1与Lagrange分析完全类似§3HermiteInterpolation一般地，已知x0

,…,xn

处有y0

,…,

yn

和y0’

,…,yn’

，求H2n+1(x)满足H2n+1(xi)=yi，H’2n+1(xi)=yi’。解：设+=ni)()()(=0iixxyixH2n+1n=0iyi′其中

i(xj)=ij,i′(xj)=0,

(xj)=0,

(xj)=ij

ii(x)有根x0

,…,xi,…,xn且都是2重根由余下条件i

(xi)=1和

i′(xi)=0可解Ai

和Bi

(x)i有根x0

,…,xn,除了xi

外都是2重根i)()(iili2(x)xxCx-=又i(xi)=1Ci

=1i)(x)(ili2(x)xx-=这样的Hermite插值唯一i′)()()(2xlBxAxiii+=in+1次所以类似有：,则证明略.§3HermiteInterpolationQuiz:给定xi=i+1,i=0,1,2,3,4,5.下面哪个是i(x)的图像？x0--10.5123456yxy0---10.5123456斜率=1

求Hermite多项式的基本步骤：写出相应于条件的i(x)、i(x)的组合形式；对每一个i(x)、i(x)找出尽可能多的条件给出的根；根据多项式的总阶数和根的个数写出表达式；根据尚未利用的条件解出表达式中的待定系数；最后完整写出H2n+1(x)。注：待定系数法仍适用，但插值节点多时比较麻烦。1分析（方法1）：误差：#例2:建立插值多项式)(3xH,使之满足插值条件

îíì¢===)()(2,1,0)()(00'33xfxHixfxHii.

方法2：（用带有重节点的差商表）#§7分段低次插值

/*piecewisepolynomialapproximation*/

为了保证插值函数的逼近效果，需要较多的插值节，导致较高的多项式次数。然而在实际应用中,很少采用高次插值:①在两相邻插值节点间,插值函数未必能够很好地近似被插值函数;②对于等距节点的牛顿插值公式,函数值的微小扰动可能引起高阶差分有很大的变化。7-1多项式插值的问题RememberwhatIhavesaid?IncreasingthedegreeofinterpolatingpolynomialwillNOTguaranteeagoodresult,sincehigh-degreepolynomialsareoscillating.例：在[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

n越大，端点附近抖动越大，称为Runge现象Ln(x)f(x)分段低次插值分别采用10次、15次、20次的等距节点插值多项式。随着插值次数的提高,在范围内的近似程度并没有变好,反而变坏.高次插值并不一定带来更好的近似效果。§4PiecewisePolynomialApproximation7-2分段线性插值

/*piecewiselinearinterpolation*/在每个区间上，用1阶多项式(直线)逼近f(x):记，易证：当时，一致失去了原函数的光滑性。

分段Hermite插值

/*Hermitepiecewisepolynomials*/给定在上利用两点的y及y’构造3次Hermite函数导数一般不易得到。Howcanwemakeasmoothinterpolationwithoutaskingtoomuchfromf?Headache…7-3分段三次Hermite插值

/*Hermitepiecewisepolynomials*/§8三次样条

/*CubicSpline*/

分段插值法具有一致的收敛性,但它只保证插值函数整体的连续性,但在连接处不一定光滑，不能够满足精密机械设计（如船体、飞机、汽车等的外形曲线设计）对函数光滑性的要求。早期的工程技术人员在绘制给定点的曲线时，使用一种具有弹性的细长木条（或金属条），称之为样条（Spline），强迫它弯曲通过已知点。弹性力学理论指出样条的挠度曲线具有二阶连续的导函数，并且在相邻给定点之间为三次多项式，即为数学上的三次样条插值曲线。定义设。三次样条函数

,且在每个上为三次多项式

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

/*cubicsplineinterpolant*/.注：三次样条与分段Hermite插值的根本区别在于S(x)自身光滑，不需要知道f的导数值（除了在2个端点可能需要）；而Hermite插值依赖于f在所有插值点的导数值。f(x)H(x)S(x)三次样条插值函数在每一个小区间上是3次的多项式,在整个插值区间上有4n个系数.且有4n-2个约束:内节点

边界节点要确定4n个系数，还需附加2个约束条件.常用的约束条件有以下三类：此时一般有成立.③.周期性边界条件,②.弯矩边界条件特别的称为自然边界条件.①.转角边界条件

构造三次样条插值函数的三转角方程

/*methodofbendingmoment*/

构造三次样条插值函数的三转角方程

/*methodofbendingmoment*/

构造三次样条插值函数的三转角方程

/*methodofbendingmoment*/

构造三次样条插值函数的三转角方程

/*methodofbendingmoment*/

构造三次样条插值函数的三转角方程

/*methodofbendingmoment*/

构造三次样条插值函数的三转角方程

/*methodofbendingmoment*/§5CubicSpline

构造三次样条插值函数的三弯矩法

/*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)(21121S[j]’(x)=jjjjjjjjBxAhxxMhxxM++-+---6)(6)(3131S[j](x)=利用已知S[j](xj1)=yj1

S[j](xj)=yj

可解§5CubicSplinejjjjjjjhMMhyyA611-----=jjjjjjjjjjjjhxxhMyhxxhMyBxA12211)6()6(-----+--=+下面解决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)=S[j+1]’(xj)，合并关于Mj1、

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

j

1n1即：有个未知数，

个方程。n1n+1还需2个边界条件

/*boundaryco