Chapter4插值与拟合-ZTY

1、1第四章 插值与拟合 在科学研究和工程中，常常会遇到计算函数y=f(x)的函数值、导数值、零点、极值或积分值等一类问题。然而函数关系往往是很复杂的，甚至没有明显的解析表达式。 24.1 引例及问题综述4.1.1 问题综述3xx0 x1xnf(x)f(x1)f(x1)f(xn)例如，根据观测或实验得到一系列的数据，确定了与自变量的某些点相应的函数值。 当函数比较复杂或根本无法写出解析式时，往往根据观测数据构造一个适当的简单的函数近似地代替要寻求的函数。 当近似的指标含义不同时，就构成了插值插值与拟合拟合两种研究方法。当精确函数当精确函数 y = f(x) 非常复杂或未知时，在一系列节非常复杂或未

2、知时，在一系列节点点 x0 xn 处测得函数值处测得函数值 y0 = f(x0), yn = f(xn)，由，由此构造一个此构造一个简单易算的近似函数简单易算的近似函数 P(x) f(x)，满足条，满足条件件P(xi) = f(xi) (i = 0, n)。这里的。这里的 P(x) 称为称为f(x) 的的插值函数。最常用的插值函数是插值函数。最常用的插值函数是 多项式多项式x0 x1x2x3x4xP(x) f(x)思路思路/* Lagrange Polynomial */如何求多项如何求多项式？式？5 定义定义4.1 设函数 y = f (x) 在区间a, b上有定义，且已知f (x) 在 a

3、, b上n+1个互异节点：x0 ，x1，xn 上的函数值f (x0 )，f (x1 )，f (xn ) 。若存在一个f (x)的近似函数P(x)，满足( )( ) (0,1,2, ) (4.8)iiP xf xin插值函数插值函数被插函数被插函数(原函数）原函数）插值节点插值节点插值条件插值条件插值区间插值区间误差函数误差函数( )( )( )R xf xP x插值余项插值余项6 通常满足同一个插值条件的插值函数有许多类型，如多项式多项式、三角函数类型、指数函数类型等等。 由于多项式或分段多项式计算简单，分析容易，而往往被选为插值函数。 当插值函数是多项式时称为代数插值代数插值（或多项式插多项

4、式插值值）。 即若存在一个次数不超过n次的多项式Pn(x)，满足( )( ) (0,1,2, ) (4.9)niiP xf xin则称P Pn n( (x x) )为f (x) 的n n次插值多项式次插值多项式。本章只讨论代数插值。满足条件：满足条件：P(xi)= f(xi)=yixi xj互异的点互异的点对于一般的对于一般的n n次代数插值多项式可以写成下列形式：次代数插值多项式可以写成下列形式： P Pn n( (x x)=)=a a0 0+a+a1 1x+ax+a2 2x x2 2+ + +a an nx xn n (1) (1) 使其在给定的使其在给定的n+1n+1个个互异互异的插值基

5、点上满足插值原则的插值基点上满足插值原则 P Pn n( (x xi i)=)=y yi i, i=0,1, i=0,1,n ,n (2) (2) 根据插值原则式根据插值原则式(2),代数多项式代数多项式(1)中的各个系数中的各个系数a0,a1,an应满足下列应满足下列n+1阶线性方程组阶线性方程组2001020002101 121112012()()()nnnnnnnnnnnnnnP xaa xa xa xyP xaa xa xa xyP xaa xa xa xy(3)9n次插值多项式Pn(x)是否存在？如何求解？插值误差或余项又如何估计？主要内容线性插值线性插值n = 1已知已知 x0 ,

6、 x1 ; y0 , y1 ，求，求xaaxP101)( 使得使得111001)(,)(yxPyxP 可见可见 P1(x) 是过是过 ( x0 , y0 ) 和和 ( x1, y1 ) 两点的直线。两点的直线。1x线性插值示意图P1(x)f(x)x0 x1y0y14.2.1 线性插值与抛物插值4.2 拉格朗日插值?线性插值线性插值n = 1已知已知 x0 , x1 ; y0 , y1 ，求，求xaaxP101)( 使得使得111001)(,)(yxPyxP 可见可见 P1(x) 是过是过 ( x0 , y0 ) 和和 ( x1, y1 ) 两点的直线。两点的直线。)()(0010101xxxx

7、yyyxP 101xxxx 010 xxxx = y0 + y1l0(x)l1(x) 10)(iiiyxl00011011()1;()0()0; ()1lxlxl xl x在节点上满足在节点上满足线性插值线性插值n = 1已知已知 x0 , x1 ; y0 , y1 ，求，求xaaxP101)( 使得使得111001)(,)(yxPyxP 可见可见 P1(x) 是过是过 ( x0 , y0 ) 和和 ( x1, y1 ) 两点的直线。两点的直线。)()(0010101xxxxyyyxP 101xxxx 010 xxxx = y0 + y1l0(x)l1(x) 10)(iiiyxl0001101

8、1()1;()0()0; ()1lxlxl xl x在节点上满足在节点上满足称为拉氏基函数（线性插值基函数）称为拉氏基函数（线性插值基函数） /* Lagrange Basis */，满足条件满足条件 li(xj)= ij /* Kronecker Delta */当当i=j时时 ij=1；当当ij时时 ij=0 抛物线插值抛物线插值n = 2可见可见 P2(x) 是过是过 ( x0 , y0 ) 、 ( x1, y1 )、 ( x2, y2 )的抛物线的抛物线 。使得使得已知已知 x0 , x1 , x2 ; y0 , y1 ，y2 ，求，求22012()Pxaa xa x111001)(,

9、)(yxPyxP 222)(yxP 000102101112202122()1;()0;()0()0; ()1; ()0()0;()0;()1lxlxlxl xl xl xlxlxlx 在节点上满足在节点上满足基函数都是二次的函数基函数都是二次的函数2( )( )iiyP xl x y基函数方法基函数方法为了满足插值条件在节点上有：为了满足插值条件在节点上有：012102201( )()()( )()()()()()lxA xxxxl xB xxxxlXC xxxx每个基函数每个基函数有两个零点有两个零点?000102101112202122()1;()0;()0()0; ()1; ()0()

10、0;()0;()1lxlxlxl xl xl xlxlxlx 在节点上条件在节点上条件求系数求系数2( )( )iiyP xl x y基函数方法基函数方法1201201020211012101201220212021()()1( )()()()()()()1( )()()()()1()()( )()()()()xxxxlxAxxxxxxxxxxxxBl xxxxxxxxxxxxxClxxxxxxxxxx0 x0012102201( )()()( )()()()()()lxA xxxxl xB xxxxlXC xxxxn 1希望找到希望找到li(x)，i = 0, , n 使得使得 li(xj)

11、= ij ；然后令；然后令 niiinyxlxP0)()(，则显然有，则显然有Pn(xi) = yi 。li(x)每个每个 li 有有 n 个根个根 x0 xi xn njj i jiniiixxCxxxxxxCxl00)().().()( j i jiiiixxCxl)(11)( njijjijixxxxxl0)()()( niiinyxlxL0)()(Lagrange Polynomial164.2.2 拉格朗日插值多项式由一次、二次插值多项式的表示方法很容易推广到一般情形。如何构造通过n+1个节点x0 x1 xn的n次插值多项式Ln(x)，满足() (0,1,2, ) (4.15)njj

12、Lxyjn1 () ( ,0,1,2, ) (4.16)0 jkkjlxj knkj定义定义4.2 若n次多项式lj (x) (j=0,1,2, n)在n+1个节点x0 x1 xn上满足条件则称这n +1个n次多项式lj (x) (j=0,1,2, n)为节点x0 x10),并注意到xi=x0+ih(i=0,1,n),则可简化为 20000010000000(1)(1)(1)()2!-() ,!nnkknkjt tt tt nN xthft fffnfx xtjtffkh 其中 这个用向前差分表示的插值多项式，称为牛顿向前插值公牛顿向前插值公式。式。200000120011( )()()()2

13、! ()()()!nnnnffNxfxxxxxxhhfxxxxxxn h64由插值余项公式，可得 (1)1( )( )( )(1)!nnnfRxxn1010( )()()()()nnniixxxxxxxxx),( )()!1()() 1()(0) 1(10nnnnxxfhnntttthxR65例例 从给定的正弦函数表出发计算sin(0.12)，并估计截断误差。 取x0=0.1,此时 ，构造差分表。表中带下划线各数依次为sinx在x0=0.1处的函数值和各阶差分。2 . 01 . 01 . 012. 00hxxtxsinx0.10.20.30.40.50.60.099830.198670.295

14、520.389420.479430.56464-0.00096-0.00094-0.00091-0.00199-0.00295-0.00389-0.004800.098840.096850.093900.090010.08521f2f3fX X和和Y Y全抄下，分全抄下，分子子前面下减上，分母就是前面下减上，分母就是1 1！66 若用线性插值求sin(0.12)的近似值，则可得sin(0.12)N1(0.12)=0.09983+0.20.098840.11960用二次插值得 )00199. 0(2) 12 . 0(2 . 009884. 02 . 009983. 0)12. 0()12. 0s

15、in(2 N =N1(0.12)+0.00016=0.11976 用三次插值得 sin(0.12)N3(0.12) 11971. 0)00096. 0(6)22 . 0() 12 . 0(2 . 0)12. 0(2 N200000(1)(1)(1)()2!nnt tt tt nN xthft fffn 67 因N3(0.12)与N2(0.12)很接近，且由差分表可以看出，三阶差分接近于常数（即4f0接近于零），故取N3(0.12)=0.11971作为sin(0.12)的近似值，此时由余项公式可知其截断误差000002. 0 )4 . 0sin() 1 . 0(24)32 . 0()22 . 0

16、() 12 . 0(2 . 0)12. 0(43R),( )()!1()() 1()(0) 1(10nnnnxxfhnntttthxR68n次牛顿（牛顿（Newton）插值多项式插值多项式 差商表：差商表：牛顿向前牛顿向前插值公式插值公式 差分表：差分表：00100120101011( )(),(),()() ,()()()nnnNxf xf x xxxf x x xxxxxf x xxxxxxxxX X和和Y Y全抄下，分全抄下，分子子前面下减上，分母前面下减上，分母X X大减小！大减小！20000010000000(1)(1)(1)()2!-() ,!nnkknkjt tt tt nN x

17、thft fffnfx xtjtffkh 其中X X和和Y Y全抄下，分全抄下，分子子前面下减上，分母就是前面下减上，分母就是1 1！69 4.5 分段线性插值分段线性插值704.5.1 高次插值的病态分析 适当地提高插值多项式的次数，有可能有可能提高计算结果的准确程度。 但是决不可由此得出结论，认为插值多项式的次数越高越好，利用被插函数节点信息越多，误差越小。由插值多项式的截断误差公式可见： (1)1( )( )( )( )( )(1)!nnnnfRxf xpxxn其中其中 1010( )()()()()nnniixxxxxxxxx 截断误差同f(n+1)()与n+1(x)有关，其绝对值不一

18、定随次数n增加而减小。71龙格现象龙格现象考虑一个典型的例子，设21( ) (55)1f xxx 取等距节点105 (0,1, )kinxkn 所构造的拉格朗日插值多项式为2001()1nninjijjiijxxLxxxx在-5，5上取一点55nx 计算出n=2,4,20的( )( )nL xR x和 ninijjijijniniiiiiiniiiniiinyxxxxxxxxxxxxxxxxxxxxyyxlxL0001101100)( )()()()()()()()(72n20.137931 0.759615-0.62168440.066390 -0.3568260.42321660.0544

19、63 0.607879-0.55341680.049651 -0.8310170.880668100.047059 1.578721-1.531662120.045440 -2.755002.800440140.044334 50.332742-5.288409160.043530 -10.17386710.217397180.042920 20.123671-20.080751200.042440 -39.95244939.994889( )nL x ( )R x ( )f x 随着n的增加， 不但没有减少，反而成倍地增加。( )nR x 这个例子最早是在20世纪初由龙格(Runge)研究的

20、。这种节点变密但误差增大的现象称为龙格(Runge)现现象。象。734.5.2 分段线性插值 当插值节点很多时，使用高次插值多项式未必能够得到好当插值节点很多时，使用高次插值多项式未必能够得到好的效果。的效果。 我们介绍分段插值法，就是将插值区间分成若干个子区间，然后在每个子区间上使用线性插值多项式，分段线性插值分段线性插值: :通过插值点用折线段连接起来逼近通过插值点用折线段连接起来逼近f (x) 。 设f (x)在n+1个节点 a = x0 x1 x2 xn= b上的函数值为 f (xk) ( k = 0,1,2, n)，将a,b分成n个小区间，在每个小区间xk , xk+1上作线性插值，

21、得74111111( )( )() ()kkkkkkkkkkx xx xL xf xf xxxxxxxx 若要计算点xxi处函数值f (x)的近似值。可先选取两个节点xi-1与xi，使xxi-1,xi；然后在小区间上作线性插值，计算函数值。111111( )( )()( ), iiiiiiiiiix xx xf xL xf xf xxxxxxxx ninijjijijniniiiiiiniiiniiinyxxxxxxxxxxxxxxxxxxxxyyxlxL0001101100)( )()()()()()()()(75在几何上就是用通过曲线n+1个点(xi, yi)的折线，去近似代替曲线 。故分

22、段线性插值又称折线插值折线插值。76分段线性插值的余项估计式：11112( )( )()()() (,)kkkkkkf xL xfxxxxx x注意在xi-1, xi上21114()(), kkkkkkxxxxhhxx推得)()()(10211xxxxfxR 1111112218max( )( )max()()() max( )kkkkkkkkkxx xxx xkxx xf xL xfxxxxhfx X轴之下的开口向上的抛物线77这表明，只要小区间长度hk充分小，便可保证 充分靠近f (x) ，即 时分段线性插值函数 收敛于被插函数f (x) 。1( )Lx1max0kknhh1( )Lx01

23、10110121801max( )( )max( )( ) max max( )( ) max max( )( ) maxmnkkkka x bxx xk nxx xk nxx xkk nf xL xf xL xf xL xf xL xh 121801ax( ) max( ) (max)kkxx xka x bk nfxhfxhh 可见，分段线性插值的余项只依赖于二阶导数的界。78 在科学实验和生产实践中，经常要从一组实验数据（xi,yi）(i=1,2,m)出发，寻求函数y = f (x)的一个近似表达式y =p(x)（称为经验公式）。从几何上看，就是希望根据给定的m个点（xi,yi） ，求曲

24、线y = f (x)的一条近似曲线y =p(x) 。因此，这是一个曲线拟合的问题。 多项式插值虽然在一定程度上解决了由函数表求函数的近似表达式问题，但用它来解决这里提出的问题，是有明显缺陷的。 4 4.7 .7 曲线拟合的曲线拟合的最小二乘法最小二乘法79多项式插值的缺陷多项式插值的缺陷 由实验提供的数据通常带有测试误差。如果要求近似曲线y=p(x)严格地通过所给的每个数据点（xi,yi），就会使曲线保留着原有的测试误差。当个别数据的误差较大时，插值效果显然是不理想的。 怎样从给定的一组实验数据出发，在某个函数类中寻求一个“最好”的函数p(x)来拟合这组数据，是一个值得讨论的课题。 随着拟合效

25、果“好”、“坏”标准的不同，解决此类问题的方法也就不同。这里，将介绍一种最常用的曲线拟合方法，即最小二乘法。 80(a)(b)(c)(d)例例81多项式插值：必须通过所有的数据点。 最小二乘法：不需要通过所有的数据点，而是数据点误差之和最小。82什么是最小二乘法什么是最小二乘法 在一般情况下，不能要求近似曲线y=p(x)严格地通过所有数据点（xi,yi）,亦即不能要求拟合函数在xi处的残差（亦称偏差） ), 2 , 1( )(miyxpiii都严格地等于零。但是，为了使近似曲线能尽量反映所给数据点的变化趋势，要求 都较小还是需要的。达到这一目标的途径很多，常见的有： i1）选取p(x),使残差

26、绝对值之和最小，即min)(11miiimiiyxp83（2）选取p(x),使残差最大绝对值最小，即 min)(maxmax11iimimiiyxp（3）选取）选取p(x),使使残残差平方和差平方和s最小最小，即 min)(2112iimimiiyxps 为了便于计算、分析与应用，较多地根据“使残差平方和使残差平方和最小最小”的原则（称为最小二乘原则）的原则（称为最小二乘原则）来选取拟合曲线y=p（x）。 按最小二乘原则选择拟合曲线的方法，称为最小二乘法最小二乘法（Least-Squares MethodLeast-Squares Method）。 p（x）称为m个数据（xi,yi）(i=1,

27、2,m)的最小最小二乘法拟合函数二乘法拟合函数。y p(x)近似反映了变量x和y之间的函数关系 y = f (x)，称为经验经验公式公式或数学模型数学模型。 f (x)称为被拟合函数被拟合函数。84定义4.7 设x1, x2,xn为互不相同的点，01( ),( ),( )1mxxxm是个已知函数。如果存在不全为零的常数 ，使得 01,mC CC0 01 1( )( )( ) 0 (1,2, ),jjm mjCxCxCxjn则称函数 （关于节点x1, x2,xn ）是线线性相关性相关的，否则称为线性无关线性无关。01( ),( ),( )mxxx简单的说，如果线性相关，则其中任一函数都可以用其他

28、函数的线性关系来表示85定义定义4.8 4.8 给定数据(xj ,yj) (j=1,2,n)。xx1x2xny=f (x)y1y2yn( )(0,)kx km其中2201110(,)()()(4.63)nnmmjjkkjjjjkaaaP xyaxy 0011( )( )( )( ) (4.62)mmP xaxaxax设拟合函数的形式为为已知的线性无关函数，求系数01,ma aa使得最小，若22*10100()min() (4.64)knmnmkkjjkkjjaRjkjkkmaxyaxy 86*0011( )( )( )( )mmP xaxaxax则称相应的为最小二乘拟合函数最小二乘拟合函数特别

29、地，若*2*012( )mmP xaa xa xa x则称P(x)为m次最小二乘拟合多项式次最小二乘拟合多项式。87 用最小二乘法解决实际问题包含两个基本环节：（1） 先根据所给数据点的变化趋势与问题实际背景确定线性无关的函数 ，即确定P(x)所具有的形式。（2 ）然后按最小二乘原则(残差平方和S最小) 求取最小二乘解 ，即确定其系数 。 *( )P x*(0,1, )ka km则可以得到拟合曲线 ( )(0,)kx km*0011*0( )( )( )( )( )mmmkkkP xaxaxaxax2211( )minnhjjjjjSP xyn确定确定k(x)的形式的形式确定确定ak*88最小

30、二乘解的求法最小二乘解的求法最小二乘解 应满足条件得: 点 是多元函数 的极小点。由多元函数取极值的必要条件知，*01(,)maaa满足方程组12(,)0 (0,1,)mka aakma用求多元函数极值的方法求最小点*01(,)maaa22*10100()min() (4.64)knmnmkkjjkkjjaRjkjkkmaxyaxy *01(,)maaa2201110(,)()()(4.63)nnmmjjkkjjjjkaaaP xyaxy *01(,)maaa89将式(4.63)两边对ak求偏导，并令2201110(,)()(4)( .63)nnmmjjkkjjjjkaaaP xyaxy 21

31、1*10()() 2() =2()() =0 (0,1,.,)nnjjjjjjjkkknmiijjkjjiP xP xyP xyaaaaxyxkm 即 *00111()()()()0nkjjjmmjjjxaxaxaxy亦即 *00111111( ) ( )( ) ( )( )( )( )nnnnkjjkjjmkjmjkjjjjjjaxxaxxaxxx y90*00111111( ) ( )( ) ( )( )( )( )nnnnkjjkjjmkjmjkjjjjjjaxxaxxaxxx y*011()()() (0,1,) (4.65)mnnijkjijkjijjxxayxkm 即为了将上式表达

32、得更简洁，引进内积记号内积记号，在线性代数中，Rn中两个向量 的内积定义为11( ,)(,)TTnnxxyyxy和11( , )nnx yx yx y将它加以推广有下面的定义。91定义定义4.9 设u(x)与v(x)是两个已知函数，记1212( (), (), () ,( (), (), () ,TTnnu xu xu xv xv xv xuv令1( , )() () (4.66)njjju x v xu v称之为u和和v的内积的内积。内积的性质：121) ( , )0,( , )00;2) ( , )( , );3) (,)( ,)( ,);( (), (), ()4) (, )( , )(

33、Tnw xw xw xu uu uuu vu vuv wu wv wwu vu v当且仅当其中为任意实数）(u,v) = (v,u);92*0011(,)(,)(,)( ,) (0,1, )kkmkmkaaaykm 利用内积的定义，则式*00111111( ) ( )( ) ( )( )( )( )nnnnkjjkjjmkjmjkjjjjjjaxxaxxaxxx y其中01020010()()()()nnxxxx11121111()()()()nnxxxx112211()(),()()mmmmnnmnnxyxyyxyxy93写成矩阵形式即(？) *0010000*0111111*01( ,)(

34、 ,)(,)( ,)( ,)( ,)(,)( ,) (4.67)( ,) ( ,)(,)( ,)mmmmmmmmyayaya 法方程组法方程组 或正规方程组或正规方程组, ,其中系数矩阵是对称的。其中系数矩阵是对称的。*0011(,)(,)(,)( ,) (0,1, )kkmkmkaaaykm 当函数 线性无关时，方程组是对称正定的，因此有唯一解。01( ),( ),( )mxxx求出(4.67)的解后，代入(4.62)即可得最小二乘拟合函数。0011( )( )( )( ) (4.62)mmP xaxaxax94代数多项式拟合am9501( ) mmP xaa xa x11102111111

35、121111nnnmjjjjjjnnnnmjjjjjjjjjnnnnnmmmmjjjjjjjjjnxxyaxxxx yaaxxxx yam96例例 1 1 某种铝合金的含铝量为 x%，其熔解温度为 y， 由实验测得 x 与 y 的数据如表所示。试用最小二乘法 建立 x 与 y 之间的经验公式。 i xi yi xi2 xiyi 1 2 3 4 5 6 36.9 46.7 63.7 77.8 84.0 87.5 181 197 235 270 283 292 1 13 36 61 1. .6 61 1 2 21 18 80 0. .8 89 9 4 40 05 57 7. .6 69 9 6 60 05 52 2. .8 84 4 7 70 05 56 6. .0 00 0 7 76 65 56 6. .2 25 5 6 66 67 78 8. .9 9 9 91 19 99 9. .9 9 1 14 49 96 69 9. .5 5 2 21 10 00 06 6. .0 0 2 23 37 77 72 2. .0 0 2 25 55 55 50 0. .0 0 3 39 96 6. .6 6 1 14 45 58 8 2 28