第三章插值与拟合

1、L/O/G/O第三章第三章 插值与拟合插值与拟合 在科学研究与工程技术中，常常需要从一组测量数据出发，寻找变量之间函数关系的近似表达式，这就是函数逼近问题，插值和拟合都是函数逼近的重要方法。 本章主要讲解关于插值和拟合的相关知识，包括对应的原理，计算公式的构造和误差估计等等。 应用举例:胡克定律-弹簧在力F的作用下伸长x，在一定范围内服从胡克定律，F与x成正比，即F=kx，k为倔强系数，现在实验中测得以下一组x，F的数据，并在(x，F)坐标下作图。可以看出，当F达到一定数值后，就不服从胡克定律。 试由数据确定 k，并给出不服从胡克定律时的近似公式。1 .216 .206 .198 .186 .

2、157 .116 . 69 . 35 . 1)(1715131297421)(kgFcmx0246810121416180510152025F (kg)x (cm)观察一下数据分布图，该分布满足什么趋势？第一节第一节 插值问题插值问题 在生产实践的许多领域里，例如机械工业、造船、汽车制造，常常有这样的问题：给了一批离散样点，要求作出一条光滑曲线（乃至于曲面），使其通过或尽可能的靠近这些样点，以满足设计要求或者据此进行机械加工这就是插值问题。 插值法历史悠久。据考证在公元六世纪时，我国的刘焯（zhuo）已把等距二次插值应用于天文计算。十七世纪时，Newton和Gregory（格雷格里）建立的等距

3、节点上的一般插值公式，十八世纪时，Lagrange给出了更一般的非等距节点插值公式。 刘焯（544-610）隋天文学家。字士元。着力研习九章算术、周髀、七曜历书等；著有稽极10卷，历书10卷。提出新法，编有皇极历，在历法中首次考虑太阳视差运动的不均匀性，创立用三次差内插法来计算日月视差运动速度，推算出五星位置和日、月食的起运时刻，这是中国历法史上的重大突破。插值问题的基本概念 设函数y=f (x)在区间a,b上有定义， 且它在该区间上的n+1个互异点 的函数值已知，即 如果选取简单函数P(x)作为y=f (x)的 近似表达式，并满足以下条件： （3-1）bxxxxan210iiyxf)(ni,

4、 1 , 0 iiyxP)(ni, 1 , 0 这样的函数近似问题就称为插值问题。 称为插值条件；满足插值条件的近似函数P(x)就称为f (x)的插值函数，而y=f (x)称为被插值函数；互异点x0,x1,.,xn被称为插值节点(简称节点)，x称为插值点；区间a,b称为插值区间。 这样对f (x)在区间a,b的各种运算就用对插值函数P(x)的运算取而代之。iiyxP)(ni, 1 , 0 代数插值的几何意义，即通过n+1个点 做代数曲线近似代替原曲线y=f (x)。), 2 , 1 , 0)(,(niyxii)(xfyO)(xPy xx1ynxny1插值函数是否存在？插值函数是否唯一？如何表示

5、插值函数？如何估计被插值函数与插值函数间的误差？构造插值函数需要关心下列问题：插值多项式存在的唯一性定理2-1 在n+1个互异节点(xi，f (xi)上满足插值条件 的次数不高于n次的插值多项式 nnnxaxaxaaxP2210)(iiyxP)(ni, 1 , 0 （3-2）存在且唯一。)()()(22101121211000202010nnnnnnnnnnxfxaxaxaaxfxaxaxaaxfxaxaxaa未知量a0 , a1 , .an的系数行列式称为范德蒙（Vandermonde）行列式（3-3）证明：根据插值条件,（3-2）式中的系数a0 , a1 , .an应满足以下n+1阶线性方

6、程组 由于节点互异，即 ，所以 ，由克莱姆法则可知方程组（3-3）有唯一的一组解a0 , a1 , .an ，也就是插值多项式（3-2）存在且唯一。0V)(jixxji0212110200)(111jinjinnnnnnxxxxxxxxxxxV（3-4）第二节第二节 拉格朗日插值多项式拉格朗日插值多项式 拉格朗日插值多项式是一种在形式上不同于（3-2）的插值多项式，这种多项式的明显优点是无需求解方程组，只要给出了n+1个互异节点及对应的函数值，便能直接写出这种形式的插值多项式。 拉格朗日插值多项式的构造方法 已知函数y =f (x)在 n+1个节点处的函数值为 ，现要作一个n次插值多项式 ，并

7、使在节点xi处满足nxxx10), 2 , 1 , 0)(nixfyii)(xLn)(xLn), 2 , 1 , 0()(niyxLiin将构成n+2个方程的联立方程组，而未知量只有n+1个，根据非齐次线性方程组有解的充要条件必有)(2210 xLxaxaxaannn0111)(1211211002002nnnnnnnnnyxxxyxxxyxxxxLxxx在方程组（3-3）上加上一个方程 利用范德蒙行列式的计算公式，经整理后可得到的表达形式如下： 其中 叫做拉格朗日插值多项式的插值基函数，显然它满足njiijijxlji, 2 , 1 , 0,;10)()(xli)()()()()()()(1

8、1101110niiiiiiiniiixxxxxxxxxxxxxxxxxxxxxl)()()()(1100 xlyxlyxlyxLnnn（3-5）（3-6）两个节点带公式得线性插值也可以直接列方程组求解，即)(,(),(,(1100 xfxxfx)()()()()()(101001011xfxxxxxfxxxxxLxaaxL101)()(xfy)(1xL（3-7）线性插值多项式 二次插值多项式)()()()(02010212xfxxxxxxxxxL)()()()()()(21202101012102xfxxxxxxxxxfxxxxxxxx（3-8）)(xfy)(2xL三个节点)(,(),(,(

9、),(,(221100 xfxxfxxfx插值多项式的余项（截断误差）用n次插值多项式进行计算所产生的误差叫做插值多项式的余项或截断误差，记作Rn （3-9）)()()(xPxfxRnn定理2-2设函数y=f (x)在插值区间上具有直到n+1阶的连续导数，x0，x1，.，xn是互异节点，则插值多项式(3-2)的余项为插值余项的估计： IxnfxxnfxRnnniinn),()!1()()()!1()()(1) 1(0) 1(niiniinnxxnMxxnfxR00) 1()!1()()!1()()()() 1(),(maxxfMnbax式中（3-10）（3-11） (3-10)式称为多项式插值

10、的余项公式，由于无法求出公式中 的准确值，所以实际计算中用它来估计误差仍有困难。因此，只在理论分析中应用。实际应用中通常采用(3-11)式来作为多项式插值的余项 估计式。 试用线性插值和抛物线插值两种方法分别计算ln11.75，并估计截断误差。xi1011121314yi2.30262.39792.48492.56492.6391例1. 已知y=lnx线性插值：选11,12两点带公式计算4849. 2)1112()11(3979. 2)1211()12()(1xxxL4409. 1087. 0 x4632. 2)75.11(1L余项0007. 0)75.11(75.11ln1Lx0008. 0

11、)1275.11)(1175.11(! 2111)(21xRniiniinnxxnMxxnfxR00) 1()!1()()!1()()(0006. 0)(2xRniiniinnxxnMxxnfxR00) 1()!1()()!1()()(5649. 2)1213)(1113()12)(11(4849. 2)1112)(1312()11)(13(3979. 2)1311)(1211()13)(12()(2xxxxxxxL二次插值：选11,12,13三点带公式计算4638. 2)75.11(2L余项0001. 0 x例2:已知e-x在x=1,2,3的值如下，试用二次Lagrange插值求e-2.1的

12、近似值，并进行误差估计。x123e-x0.3678794410.1353352830.049787068解：二次插值多项式为3212) 23)(13 () 2)(1() 32)(12 () 3)(1() 31)(21 () 3)(2()(exxexxexxxL因 ，故误差的余项2121 . 2)32)(12()31 . 2)(11 . 2()31)(21 ()31 . 2)(21 . 2() 1 . 2(eeLe1 21)(max eexx00607001. 0)31 . 2)(21 . 2)(11 . 2(! 3) 1 . 2(12eR120165644. 0)23)(13()21 . 2)

13、(11 . 2(3eniiniinnxxnMxxnfxR00) 1()!1()()!1()()(第三节牛顿插值多项式第三节牛顿插值多项式 当使用拉格朗日插值多项式为了提高精度增加插值节点时，插值多项式就得重新构造，整个公式改变之后，之前的计算结果不能继续使用，计算工作必须全部从头做起。为了克服这一缺点，下面介绍另一种插值多项式牛顿插值多项式，它的使用比较灵活，当增加插值节点时，只要在原来的基础上增加部分计算工作量，而原来的计算结果仍可利用，这为实际计算带来了方便。差商的定义定义 已知函数f (x)在n+1个互异节点xi (i=0,1,2,.n)上的函数值分别为f (xi)(i=0,1,2,.n

14、)称为f (x)关于节点 的一阶差商。称为f (x)关于节点 的二阶差商。iiiiiixxxfxfxxf111)()(,1,iixxiiiiiiiiixxxxfxxfxxxf212121,21,iiixxx则把一般地，称为f (x)关于节点 的k阶差商。当k=0时称为 关于节点 的零阶差商，记为 。ikikiiikiiikiiixxxxxfxxxfxxxf11211,kiiixxx,1)(ixfix ixf（3-12）性质1 函数f (x)关于节点x0,x1,.,xk的k阶差商f x0,x1,.,xk可以表示为函数值f (x0), f (x1),.,f (xk)的线性组合，即kiikikxxf

15、xxxf0110)()(,例3：021021210,xxxxfxxfxxxf)()()()()()(120222101120100 xxxxxfxxxxxfxxxxxf差商的性质kxxxf,10kjkjjjjjjjjxxxxxxxxxxxf01110)()()()(分析 ：用归纳法证明证明：（1）当k =1时, 010110)()(,xxxfxfxxf 011100)()(xxxfxxxf 时成立，即有假设当1 nknjnjjjjjjjnxxxxxxxxxfxxxf111121)()()()(，101110110)()()()(njnjjjjjjjnxxxxxxxxxfxxxf，011021,

16、xxxxxfxxxfnnn111110)()()(1)(njnjjjjjjnjxxxxxxxxxxxf#)()()()(1011110njnjjjjjjjjxxxxxxxxxxxf)()()()()()()(110201001100nnnnnnnxxxxxxxfxxxxxxxfxx)()()(1)()(11100njjjjjjnjxxxxxxxxxxxf,10nxxxf由差商定义性质3 设f (x)在包含互异节点x0,x1,.,xk的闭区间a,b上有n阶导数，则n阶差商与n阶导数之间有如下关系banfxxxfnn,!)(,)(10（3-13）性质2 差商与其所含节点的排列次序无关。iiiiii

17、iiiiiiixxxfxxxfxxxfxxfxxf,12212111利用差商的递推定义，差商可以用列差商表的方法来计算：造差商表（差商的计算）33221100fxfxfxfx3,22,11,0 xxfxxfxxff x0,x1,x2,x33,2,12,1,0 xxxfxxxf牛顿插值多项式及其余项由差商定义得出满足插值条件Nn(x)=f (xi)(i=0,1,2,n)的插值多项式Nn(x)。设xi (i=0,1,2,.n)为n+1个插值节点，且 ，则由差商定义有bax,), 2 , 1 , 0(nixxi000)()(,xxxfxfxxf移项整理可得再利用二阶差商的定义和性质2可得移项得将此式

18、带入（3-14）,得)(,)()(000 xxxxfxfxf110010,xxxxfxxfxxxf)(,110100 xxxxxfxxfxxf)(,)(,)()(10100100 xxxxxxxfxxxxfxfxf（3-14）重复以上过程可得记即n次牛顿插值多项式，由得其对应余项为)()(,)()(,)(,)()(101011010102100100nnnnxxxxxxxxxxfxxxxxxxxxfxxxxxxxfxxxxfxfxf)()(,)(,)()(110100100nnnxxxxxxxxxfxxxxfxfxN)()(,)(010nnnxxxxxxxxfxR)()()(xRxNxfnn（

19、3-15）由插值多项式的存在和唯一性定理，可得所以当f (x)在(a,b)上有n+1阶导数时，两者的余项也相等，即由此可得这就是差商的性质3。)()!1()()(,)(1) 1(10 xnfxxxxfxRnnnnn),(,)!1()(,) 1(0banfxxxfnn)()(xLxNnn（3-16）例4：已知一组观察数据为试用此组数据构造3次牛顿插值多项式N3(x)，并计算N3(1.5)的值。解：先按数据表造出差商表i0123xi1234yi0-5-6343213650915521相应的牛顿插值多项式只需将表中加下划线的数字代入公式即可或整理得) 2)(1( 2) 1( 50)(3xxxxN)

20、3)(2)(1(xxx34)(233xxxN625. 235 . 145 . 1) 5 . 1 (233N) 25 . 1)(15 . 1 ( 2) 15 . 1 ( 50) 5 . 1 (3N) 35 . 1)(25 . 1)(15 . 1 (625. 2例5：已知f (x)=sh(x)数值表如下：用Newton插值公式求f (0.596)的近似值。xk0.400.550.650.800.90f(xk)0.410750.578150.696750.888111.02652解：先造差商表1.38411.18601.27571.11600.43360.35880.28000.2140.1970.

21、0340.800.900.650.550.400.888111.026520.696750.578150.41075将x=0.596代入得：由牛顿公式得四次插值多项式为：N4=0.41075+1.1160(x-0.40)+0.2800(x-0.40)(x-0.55)+0.197(x-0.40)(x-0.55)(x-0.65)+0.034(x-0.40)(x-0.55)(x-0.65)(x-0.80)63192. 0)596. 0()596. 0(4 Nf插值的上机算法1.Lagrange插值法(1).输入 ,令 （输入节点数值,插值式赋初值0）(2).对 计算 （计算基函数） （累加求插值式）

22、),.,(,n210iyxii0 xLn)(,.,n210i nij0jjijixxxxxl,)(iinyxlxLxL)()()(2.Newton插值法(1).输入 （输入节点数值,插值式阶数）(2).对 计算 f x0,x1,.,xk （计算各阶差商）(3).计算 Nn(x)=f (x0)+f x0,x1(x-x0)+. +f x0,x1,.,xn(x-x0)(x-x1).(x-xn-1) （累加求插值式）);,.,(,n210iyxniik21in321k,.,.,多项式插值拉格朗日插值牛顿插值三种插值方法的比较nnnxaxaxaaxP2210)()()()()(1100 xlyxlyxl

23、yxLnnn)()()()()()()(11101110niiiiiiiniiixxxxxxxxxxxxxxxxxxxxxl)()(,)(,)()(110100100nnnxxxxxxxxxfxxxxfxfxN多项式插值多项式插值拉格朗日插值多拉格朗日插值多项式项式牛顿插值多项牛顿插值多项式式优点方法简单，易于理解无需计算形式便于记忆增加节点时原有的项不变缺点节点很多时，计算量大节点增加时，整个公式改变先计算差商，再写公式差分和等距节点公式差分与差商的区别等距节点),.,2 , 1 , 0(0niihxxihxxii)(1iiiyyy11iiiyyy)()(12iiiiyyyyiiiiiyyy

24、yy1212步长向前差分向后差分二阶差分L/O/G/O第四节第四节 分段插值分段插值龙格现象 用插值多项式逼近被插值函数时，为了达到减小余项绝对值的目的，采用提高插值多项式最高次数的方法：例6：给定函数 在区间-1，-1上取插值节点为建立插值公式22511)(xxfninixi, 2 , 1 , 0,21niiinxlxfxL1)()()(-1010.00.51.01.52.0 XY 图上的实线表示被插值函数，虚线表示插值函数；从图看出在x=-0.96及x=0.96附近，插值函数L10(x)对于被插值函数f(x)的误差很大，这种现象称为龙格(Runge)现象。当取n=10时，将插值与被插值函数

25、作图如下： 我们从分析插值余项的变化来分析这种现象的产生原因。余项公式：niinnxxnfxR0) 1()()!1()()( 许多函数的高阶导数的绝对值随着导数阶数的增加而迅速增加，因而余项的绝对值随之迅速增加。xxf1)(1)(1!) 1()(nnnxnxf例6： 另一点是插值节点差值的增加，对于固定的x，互异节点中只有少数与x邻近的节点与它的差的绝对值较小。随着节点的增加，与其差的绝对值也会随之增大。 综合以上两点我们发现用提高插值函数阶数来减小插值的余项的方法是不可行的。实际中使用的是一、二、三次插值。分段低次插值 人们从龙格现象中发现高次插值并不能减小误差，从而寻找能减小误差的另外的插

26、值方法，分段插值就是其中一种。 将给定的插值区间a,b做分割，在分割后的每一段小区域内使用低次插值，称为分段插值。分段线性插值分段二次插值折线代替曲线分段抛物线代替曲线分段线性插值分段线性插值的几何意义：用折线代替曲线把给定的插值区间a,b作分割a = x0 x1. xn= b，在每个小区间xj,xj+1上作f (x)以xj,xj+1为节点的线性插值，则在每一个对应的小区间xj,xj+1(h= xj+1-xj)上插值函数可以写为x0 xnxj区间a,b上的整个插值函数用折线方程可表示为njjjhyxlxI0)()(1111)(jjjjjjjjhyxxxxyxxxxxI拉式分段线性插值由线性插值

27、余项公式，当x xj,xj+1时有)(! 2)()(1)2(jjxxxxfxR所以)(2)(1jjxxxxMxR2121)(8)(412jjjjxxMxxMbxaxfM , )(max分段线性插值算法简单，只要分段区间充分小就能保证误差要求。其显著优点是它的局部性质，修改了某节点(xj,f (xj)的值后，仅相邻的两个区间xj-1,xj xj,xj+1受影响，缺点是插值节点处不光滑。 于是，当区间分割加密，h= xj+1-xj0时，分段线性插值收敛于f (x)。（3-17）例7：对下列数据作分段线性插值，并计算 f (1.2) 和 f (3.3)。xi-3-1239f (xi) 1251612

28、解:1111)(jjjjjjjjhyxxxxyxxxxxI,1jjxxx 2 , 12 . 10667. 211212 . 152122 . 1)2 . 1 (hI 9 , 3 3 . 3 3 . 612633 . 36693 . 3) 3 . 3 (hI因为同理22511)(xxf例8：设 ，-1 x 1，将区间-1,1 10等份，用分段线性插值近似计算f (-0.96)。解：(1)插值节点为xj=-1+ j/5 (j=0,1,10)，h= 1/5因为-0.96-1,-0.8，取此区间为线性插值区间，其上的插值函数为801. x所以 f (-0.96) Ih (-0.96) = 0.0425

29、31111)(jjjjjjjjhyxxxxyxxxxxI2 . 01) 8 . 0(2 . 08 . 0) 1()(xfxfxIh) 1(2941. 0)8 . 0(1923. 0 xx样条插值简介 如果既需要建立低次插值函数以保证稳定性，又要插值函数具有良好的光滑性，就需要采用样条函数（spline function），样条函数是指一类分段光滑，并且在各段的交接处也具有一定光滑性的函数。 样条最初是指工程绘图人员为了将一系列指定点连接成一条顺滑的曲线所采用的附有弹性的细木条或钢条，由这些样条构成的曲线在连接处能够具有连续的曲率。 1946年，美国学者肖恩博格（Schoenberg）首先提出了

30、样条函数这一概念。随后，样条函数在外形设计等领域取得了成功的应用，并与计算机辅助设计领域紧密结合，形成了一个新的交叉学科。 在计算机科学的计算机辅助设计和计算机图形学中，样条通常是指分段定义的多项式参数曲线。由于样条构造简单，使用方便，拟合准确，并能近似曲线拟合和交互式曲线设计中复杂的形状，样条是这些领域中曲线的常用表示方法。 第五节第五节 函数拟合的最小二乘法函数拟合的最小二乘法插值多项式的缺陷：测量数据带有误差，影响了插值多项式的逼近精度数据点增加时，计算量与得到的插值多项式的次数都增加 从给定实验数据出发，寻求已知函数的一个逼近函数，使逼近函数从总体上与已知函数的偏差按某种方法度量能达到

31、最小而又不一定过全部数据点。 函数拟合的几何意义就是求一条曲线，使数据点均在离此曲线上方或不远处。最小二乘法原理 用近似曲线 拟合数据(xi, f (xi)(i=0,1,2,.,m)，使得它在xi处的函数值与测量数据的偏差满足)(xmiiixfxR02)()(min即偏差的平方和最小，这就是最小二乘原理，这种拟合曲线的方法叫最小二乘法。（3-18）多项式最小二乘拟合 nnnxaxaxaaxP2210)( 用多项式 拟合数据(xi, f (xi) (i=0,1,2,.,m)，拟合数据与原数据的偏差表示如下：miininiimiiixfxaxaxaaxfx02221002)()()(由最小二乘原理

32、，得miiixfxR02)()(miininiixfxaxaxaa022210)( 利用微积分里面对函数求极值的方法， 使R对系数ai(i=1,2,.,n)的一阶导数为零，即), 2 , 1 , 0( 0niaRimiininiixfxaxaxaaaR02210001)(2:mimiininiixfxaxaxaa002210)()(对a0求导，得化简得多项式拟合系数的确定方法miiininiixxfxaxaxaaaR0221010)(2:mimiiininiiixxfxaxaxaxa00132210)()(继续对a3到an求导，得到n+1阶方程组同理，对a1求导mimiinininninimi

33、miiininiimimiininiixfxxaxaxaxfxxaxaxaxfxaxaxaa002110001210002210)()()(m+1个点n+1个方程注意：此处m和n不相等写成矩阵形式：iniiiinnininininiiiiniiiyxyxyaaaxxxxxxxxxxx1022113221利用矩阵的乘法规则，将上面的式子化为mnmnnnmnnmmmnnnmnnnmyyyxxxxxxxxaaaxxxxxxxxxxxxxxxxx212102101021211020021021011111111111将前面的方程组改写为nmmmnnxxxxxxxxxA212110200111naaa1

34、0yAAATT这个方程称为正规方程组，或法方程组myyyy21定义如下几个矩阵：（3-19）下面以线性拟合为例，直线方程为将R看成是a0,a1的二元函数，即xaaxP101)(mimiiiiinxfxaaxfxPaaR00210210)()()(),(0, 010aRaRmiiimiimiimiimiixfxaxaxxfaxam0102000100)()()()()() 1(由二元函数取极值的必要条件：得线性拟合系数的确定用解方程的方法计算系数a0和a1miiinxfxPQ02)()(线性拟合的均方误差 22220) 1()(iiiiiiixmxyxyxxa 2221) 1()() 1(iii

35、iiixmxyxmyxa例9：已知测量数据试用最小二乘法求经验直线。xi2468yi2112840建立法方程组yAAATT81614121A4028112y10aa402811286421111816141218161412110aa536811202020410aa即53612020812041010aaaa解得55. 65 .1210aa得出最后的拟合直线5 .1255. 6xy可以化为线性拟合的几种情况xbaxfln)(baxxf)(bxaxxf)(Xx lnbXaxf)(xbaxflnln)(lnbXCxF)(bxaxf)(1baXxF)(bxaexf)(bxaxf ln)(lnbxCxF)(最小二乘拟合方法的求解步骤1.用给定的数据作图2.根据数据图初步判断拟合函数的形式3.根据拟合函数的形式求解更一般地，对数据