第三章插值方法与数据拟合

1、精品文档你我共享第三章 插值方法与数据拟合所讨论的问题给复杂的函数 f ( x) 找一简单的函数 p( x) 如多项式、三角函数等，并让其满足一定的条件，让其近似的取代原函数 f ( x) 。或 有一数据表格，我们需要找一函数取近似的表征该表数据。1拉格朗日（ Lag r ang e）插值在函数类中多项式具有最简单的性质。p( x)a0 a1x1a2x2a3x3. an xn设 yf ( x) 在区间 a，b 连续的实函数已知在该区间上 n + 1个不同点 xi 的函数值 yif ( xi )i 1,2,., n或 有数据表有 n 1对数据xiyii 1,2,., n插值节点我们需要找一个 n

2、 次多项式p( x)a0 a1x1a2 x2a3 x3.an xn使得在这些点上函数值等于插值节点的值。yip( xi )1、线性插值已知两个点的函数值： x0y0x1y1做一线性函数使得在两个节点上函数值为节点值。y0p( x0 )y1p( x1)函数为：腹有诗书气自华精品文档你我共享p( x) l 0 (x) y0l1 ( x) y1xx1 yxx0 yx00x11x1x01xjxi基函数 li ( x) 为一次函数，且在节点上li ( xj)x jxi0几何意义：过两点做直线 。按 x 变化量平均。2、抛物线插值已知三个点的函数值： x0y0 x1y1 x2y2做二次函数使得在三个节点上

3、函数值为节点值。y0p( x0 ) y1p( x1) y2p( x2)函数为：p( x) l0 ( x) y0l1( x) y1l2 ( x) y2x x1 xx2 yxx0 x x2 yxx0 x x1 y2x0x1 x0x20x1x0x11x2x0 x2x1x21xjxi基函数 li ( x) 为二次函数，且在节点上li( xj)x jxi03、拉格朗日插值已知 n + 1 个点的函数值： x0y0x1y1 ,.,xnyn做 n 次函数使得在 n+ 1 个节点上函数值为节点值 。y0p( x0 )y1p( x1),.,ynp( xn )函数为：腹有诗书气自华精品文档你我共享p( x) l0

4、 ( x) y0 l1 ( x) y1 l2 ( x) y2 ,.,ln ( x) yn1x jxi基函数 li ( x) 为 n 次函数，且在节点上li (x j )x jxi0xx1xx2. xxnyx0x1 x0x2x0xn0xx0xx2xxny.x1x0 x1x2x1xn1.x x0 x x1x xi 1x xi 1.x xnyixix0 xix1.xi 1 xixi 1xi xnxi.xx0xx2xxn 1y.xnx0 xnx2xnnxn 1为了程序设计我们作如下推导：nxxjy0nxxjy1 .j 1 x0j 0 x1xjxjj1nxxjyin1xxjynj 0 xix j.0 x

5、nxjjji腹有诗书气自华精品文档你我共享nnxxjyxxi 0j 0iijj i程序结构是什么样子？DIM EN SION X(N ),Y(N )读入节点数据X= .T1= 0.0DO 100 I=0,NT2= 1.0DO 200 J= 0,NxxjT2= T2*xixj200CON TIN UET1= T1+ T2* yi100 CON TIN UEW RITE(* ,* )X,T1STOP EN D腹有诗书气自华精品文档你我共享拉格朗日插值的唯一性已知 n + 1 个点的函数值： x0y0x1y1 ,., xn yn做 n 次函数p( x)a0a1 x1a2 x2a3 x3.an xn使

6、得在 n+ 1 个节点上函数值为节点值。y0p( x0 )y1p( x1) ,.,ynp( xn )即：a0a1 x01a2 x02a3 x03.an x0ny0aa x1a x2a x3.anxny01 12 13 111aa x1ax2a x3.axny201 2223 2n2.aa x1a x2a x3.a xnyn01 n2n3 nnn求得系数 a0 , a1 , a2 , a3,., an 则多项式确定线性方程组，系数矩阵对应的行列式：1x01x02x03.x0n1x1x2x3.xn11111x1x2x3.xn02222.范德蒙（Van d er m o nd e ）行列式1x1x2

7、x3.xnnnnn所以方程有唯一解。所以多项式唯一。腹有诗书气自华精品文档你我共享2、牛顿插值公式差分、差商f ( x1 )f ( x0 )定义 4.1称 f ( x0 , x1 )x1为函数的 f ( x) 的一阶差商，x0f ( x0 , x1f ( x1 , x2 )f ( x0 , x1), x2 )x0为 f (x) 的二阶差商。x2一般地，f ( x0 , x1f ( x1 ,., xn )f ( x0 ,., xn 1),., xn )xn x0为 f ( x) 的 n 阶差商 。为统 一起见，补充定义 f ( x0 ) 为零阶差商。由差商定义，显然有f ( x0f ( x0 )

8、f ( x1 ), x1 )x1x1x0x0f ( x0 , x1 , x2 )f ( x1 , x2 )f ( x0 , x1 )x2x0f ( x0 )f ( x1 )f ( x2 )( x0 x1 )( x0x2 ) ( x1x0 )( x1 x2 ) ( x2x0 )( x2 x1 )如此类推，可以证明nf ( xj )f ( x0 , x1,., xn )j 0 ( xjx0 ).( xj xj 1 )( xj xj 1).( xj xn )腹有诗书气自华精品文档你我共享由此看出，差商的值与节点 x0 , x1 ,., xn 的排列次序无关，即f ( x0 , x1 ,., xn )

9、f ( x1 , x0 ,., xn ).f ( x1, x2 ,., xn , x0 ) 这种性质称为差商对称性。定义 4.2设函数 yf ( x) 在等距节点 xi x0 ih (i0,1,., n) 上的函数值为 f (xi )fi ，其中 h 为常数，称做步长。称fif i1fi向前差分f ififi 1向后差分fi f (xi h / 2) f (xi h / 2) f 1f 1中心差分ii22分别为 f ( x) 在 xi 处以 h 为步长的一阶向前差分，一阶向后差分和一阶中心差分。符号、分别称为向前差分算子，向后差分算子和中心差分算子。由一阶差分的定义出发，可定义二阶差分2 fi

10、fi 1fifi 22 fi 1fi ，2 fififi 1fi2 fi 1fi 2 ，2 fif1f1fi 1 2 fifi 1。i2i2一般地刻定义 n 阶差分为n fin 1 fi 1n 1 fin fin 1 fin 1 fi 1n fin 1 f 1n 1 f 1ii22腹有诗书气自华精品文档你我共享如果令 If ifi ， Efifi 1 ，则称 I 为不变算子， E 为移位算子。由于fifi 1fiEfiIf i( EI ) fi ，于是E I 。同理可得IE 1 ，E1/ 2E1/ 2牛顿插值公式有了差商的概念，前面介绍的线性插值与抛物插值可表示为N1 (x)f ( xk )f

11、 ( xk , xk 1 )( xxk ) ，抛物插值N2 ( x)f ( xk 1 )f ( xk 1 , xk )( xxk 1 )f ( xk 1 , xk , xk 1 )( xxk 1 )( xxk )事实上，按差商定义f ( x)f ( x0 )f ( x0 , x)( xx0 ) ，f ( x0 , x)f ( x0 , x1 )f ( x0 , x1, x2 )( xx1 ) ，f ( x0 , x1 , x)f ( x0 , x1 , x2 )f ( x0 , x1 , x2 , x)( xx2 )f ( x0 , x1 ,., xn 1, x)f ( x0 , x1 ,.,

12、 xn )f ( x0 , x1 ,., xn , x)( xxn )反复用后一个式子代入前面的式子，可得腹有诗书气自华精品文档你我共享f ( x)f ( x0 )f ( x0 , x1 )( xx0 )f ( x0 , x1 , x2 )( xx0 )( xx1 )f ( x0 , x1 ,., xn )( xx0 )( xx1 ).(xxn 1 )f ( x0 , x1 ,., xn , x)( xx0 )( xx1 ).( xxn )记Nn ( x)f ( x0 )f ( x0 , x1 )( xx0 ). f ( x0 , x1 ,., xn )( x x0 )( x x1 ).(x

13、xn 1)（ * ）Rn (x)f ( x0 , x1 ,., xn , x)( xx0 )( xx1 ).( xxn ) （ * ）于是f ( x)Nn ( x)Rn ( x)由于Rn (xi ) 0(i 0,1,2,., n) ，故必有Nn ( xi ) f (xi )(i 0,1,2,., n)所以式（ *）为 n 次插值多项式，由插值多项式的唯一性讨论可知Nn ( x)Ln ( x) ，且有f ( n 1)()f ( x0 , x1 ,., xn , x)。(n 1)!式（*）称作牛顿基本插值公式，它的各项系数就是函数的各阶差商，每增加一个插值节点，只需在原来的基础上多计算一项，这一性

14、质称作牛顿插值公式的承继性。腹有诗书气自华精品文档你我共享利用牛顿插值公式时，一般分两步：第一步：造差商表节函数一阶差商二阶差商三阶差商。点 值X0F(x 0)X1F(x 1)F(x0,x 1)X2F(x 2)F(x1,x 2)F(x 0,x 1,x2 )X3F(x 3)F(x2,x 3)F(x 1,x 2,x3 )F(x 0,x1,x2 ,x3)X4F(x 4)F(x3,x 4)F(x 2,x 3,x4 )F(x 1,x2,x3 ,x4)X5F(x 5)F(x4,x 5)F(x 3,x 4,x5 )F(x 2,x3,x4 ,x5)XnF(x n)F(xn - 1,xn )F(x n - 2,

15、xn - 1,x n) F(x n - 3 ,xn - 2,x n - 1,x n )第二步：利用Nn ( x)f ( x0 ) f ( x0 , x1 )( x x0 ). f ( x0 , x1 ,., xn )( xx0 )( xx1 ).(xxn 1)例如：节点表i012345Xi0 .01.02.03 .04.05.0F(xi)0 .02.012 .042.0116 .0282.0差商表XiF（xi)一阶二阶三阶四阶五阶10腹有诗书气自华精品文档你我共享222312104442301025116742240 .5628216646810 .1N5 ( x)02( xx0 )4(xx0

16、 )( xx1)2(xx0 )( x x1 )( xx2 )0.5( xx0 )( xx1 )( x x2 )( x x3 )0.1( x x0 )( x x1 )( x x2 )( x x3 )( x x4 )即：N5 ( x)02( x1)4( x1)(x2)2( x1)(x 2)( x3)0.5( x 1)(x2)( x 3)( x 4)0.1( x1)(x 2)( x3)( x4)( x5)3 、厄密（H erm i t e）插值简介高次插值逼近效果往往不理想，原因是高次函数变化复杂。所以节点的增加并不一定能提高节点之间多项式 pn (x) 与函数 f ( x) 的逼近程度。例如：函数

17、 f (x)1在区间 - 5 ，5 进行插值，我们取 n1个等分节点进行x21插值。会得到在节点上和函数值相等，节点之间误差很大。这种现象叫龙格（ Run ge ）现象 。腹有诗书气自华精品文档你我共享而利用分段低次多项式去逼近函数 f (x) 的效果要优于光滑的高次多项式 。所以可用分段线性插值。线性插值计算简单，随节点密度增加越逼近函数，但缺点是光滑度不好（导数不连续）分段三次厄密插值，就是要求函数值在节点上相等，函数的一阶导数也相等。设函数 f ( x) 在 n1个互异的节点 xii0,1,2,., n 上函数值和导数为：y0 , y1 , y2 ,., yny 0 , y 1 , y

18、2 ,., y n进行分段插值在区间 xk , xk 1 上作插值，使得在两端点函数值和导数值为：yk , yk 1y k , y k 1四个方程，我们构造三次多项式：腹有诗书气自华为三次多项式且满足：精品文档你我共享P ( x)aa xa x23012使得 P3 (xk ) ykP3 (xk 1 )P3 ( xk )y kP3 (xk 1 )a3x3yk 1y k 1令P3 ( x) yk ak (x) yk 1ak 1 ( x) y k 则 ak (x), ak 1 ( x), k (x), k 1 (x)ak (xk ) 1,ak 1 (xk ) 0, k ( xk ) 0, ak (x

19、k 1) 0, ak 1 ( xk 1 ) 1, k ( xk 1 )a k( xk )0,a k 1 (xk )0,k ( xk )a k( xk 1 )0,a k 1 ( xk1 )0, k ( xkk (x)y k 1k 1( x)k 1 ( xk )00,k 1 (xk 1)01,k 1 ( xk )01 )0,k 1 (xk 1 )1由此可求得最后 ak (x), ak1 ( x),k ( x), k 1 (x)( xxk 1)2(1 2x xk ) xk 1x xkk 0略去xkxk 1xk1xkakxxk 1)2(1 2xxk) xkx xk 1略去( x) (xk 1xkxkk

20、 nxk10x xk 1 , xk 1 腹有诗书气自华精品文档你我共享( xxk 1 )2 ( x x ) xk 1x xk 0略去xkxk 1kkk(x)( x xk 1 )2 ( x x ) xkx xk1k n略去xkxk 1k0x xk 1 , xk 14、 数据拟合问题的提出：有一组数据，数据很多，且往往有误差。需要有一函数反映这组数据。因为不同数据段误差不一样（可信度不同），同时数据多时函数中的参数有限，所以，并不要求函数通过每一个数据点。但函数尽可能科学的反映数据。往往要确定函数中的参数，使得某种统计偏差尽可能的小。或者说在某函数类中找一函数( x) ，使得某判断量最小或最大 。

21、m往往取偏差的平方和：2 yi(xi )2最小。i 1m yi(xi ) 2或2最小， Ei 为 yi 的误差。i 1Ei当函数不一样时，2 ,2 的数值不一样，这种以函数为自变量，以数为函数值的函数关系数学上叫做泛函。腹有诗书气自华精品文档你我共享yf ( x)yx2m( xi ) 22 yi( x)i 1一、用多项式进行最小二乘数据拟合现在考察 m 对实验数据或 yf (x) 的一个表函数 (xi , yi ) ， i1,2,3,., m 。函数类最简单的是多项式，于是我们首先取次数不超过 n 的代数多项式 Pn ( x)为H，即取HPn ( x) 。于是 当 y( x) Pn (x)时有

22、n( x)aa xa x2.a xna xi012i，i 0其中 a0 , a1 ,., an 取遍所有可能的实数就得到 Pn (x)的全体。n我们的目的就是要在 P x中找出一个函数 (x)a xin ( )i 0i使偏差 yi( xi) 的平方和m2(a)(a0 , a1,., an )yi( xi ) 达到极小 。i 1对于给定的一组数据 (xi , yi ) ，i1,2,3,., m 来说， (ai ) 以 a0 , a1 ,., an 为参数，也就是对应 Pi(xi ) 中一不同的 ( x) ， (a ) 有不同的值，于是问题归结为选取a0 , a1,., an 使 (a ) 取极小

23、。另外我们还要求 nm 。因 为当 nm 时即 n1m，n 次多项式由 n 1个系数所决定，而数据对只有 m 对，即使作出 ( xi )yi (i1,2,., m) 才 m 个条件，而 m n1 ，故 ( x) 至少还有一个自由度，它是不确定的，但此时的 ( x) 都使(a ) 为 0 ，当然就是极小了 。所以 对于 nm 的问题本身就没有意义 。腹有诗书气自华精品文档你我共享下面我们分几步来叙述使(a ) 达到极小的 a0 , a1 ,., an 的确定及最小的(a )值 (a0 , a1 ,., an ) 。第一步： a0 , a1 ,., an所满足的方程。m2n( x)aj x j(a

24、)(a , a ,a )yi( x )01ni1ij 0根据微积分学，极小的 a0 , a1,., an 应满足偏导数为零：mn2( yia j xi j) xi k0, k0,1,2,., n 。N+ 1 个方程aki 1j0即：mmna j x j x ky x k0,k0,1,2,.,ni iiii 1i1 j0mnmyi xi ka jxijk0,k0,1,2,.,ni 1j0i1若令mmy x kT ,x kS。i ikiki 1i 1于是有方程组na j Sj kTkk0,1,2,., n 。j 0这就是使取极小的 a0 ,a1,., an 应满足的方程。解线性方程组得到 ak ,

25、k0,1,2,., n方程的系数矩阵为腹有诗书气自华精品文档你我共享S0S1.SnS1S2.Sn1Gn 1SnSn 1.S2n我们来证明系数据正对应的行列式的值不为零。即 detGn 10 ，因为若 det Gn 10 ，则齐次方程组 Gn 1a0, (a ( a0 , a1, an )T ) ，就有非零解 a( a0 , a1,an )T，从而将齐次方程组第 k 个方程乘以 ak 并对所有 k从 0 到 n 作和有nnnnm0akaSjakajx k jj kik0j0k0j0i 1mnn)xi kj(akaji1k0j0mnak xi knj )(a j xii1k0j0mnak x j

26、k )2(k0k0nak xj k即0,j1,2,., mk 0意味着 n 次多项式，有不同的m 个零点，且 mn,所以 ak , k0,1,2,., n 全为零。二、 用函数进行最小二乘曲线拟合现在我们将前面用多项式进行拟合的情况推广到一般函数进行拟合。设 1 (x), 2 (x),.,n ( x) 是 n 个在 a,b 上定义的连续函数 。所谓函数组 1 ( x),2 ( x),.,n ( x) 是线性相关的，就是说存在 n 个不全为零的实数 a0 , a1,., an 使 a1 1 (x)a2 2 (x) . an n (x)0 在 a,b 上成立 。若上述等式在 a,b上除非 a0 ,

27、 a1,., an 全为零时成立，就说1 (x), 2 ( x),., n ( x) 是腹有诗书气自华精品文档你我共享线性无关的。现设 1( x),2 ( x),.,n ( x) 在 a,b 上线性无关的 。设 ax1x2.xmb,yf ( x) 在 xi 之值为 yif ( xi ) 。引入向量记号：( x ), (x),.,( x ) )T显然为 m 维列向量。12m所谓用函数 1 ( x),2 ( x),.,n ( x)去进行数据 (xi , yi ),i 1,2, . , m 的最小二乘曲线拟合，就是用线性拟合nQ(x)ck k ( x)k 1去逼近函数 y f ( x) ，使二者在各

28、点 xi 上偏差的加权平方和m2(c1 , c2 ,., cn )( xj ) f (x j )(c)Q( xj )*j 1n在由 1 ( x), 2 ( x),.,n ( x) 的一切可能的线性组合ck k (x) 所组成的函k 1数类 H 中为极小。若在* 中取(x)1, k ( x) xk1。于是使极小的问题就变成前面多项式数据拟合了。与前面要求 nm 一样，这里应要求 n 1m 或 nm ，否则若 n m ，则 n个 m 维向量 1, 2 ,., n ，必然线性相关。于是求* 的极小的问题就完全与求(a) 的极小一样了 。只不 过这里是函数1 , 2 ,.,n ，比前面复杂，所以引入了

29、新的内积写法。为书写简化，引进“内积”m( k , l )(x j ) k (x j ) l ( xj ) ，j 1腹有诗书气自华精品文档你我共享而称m( k , k )1/ 2(x j ) k 2 ( xj )1/ 2kj 1为 m 维列向量 k 的范数或长度。其中 ( x) 为 a, b 上给定的正连续函数，称为权函数。mn(c)(c1 , c2 ,., cn )( xj ) f (x j ) Q( xj )2ck k ( x)Q( x)j 1k1第一步：建立使(c) 取极小的 c 所满足的方程。与前面一样，应有：mnci2( xj ) y jckk (x jj 1k1mn( xj ) y

30、 jckk (x j ) i ( x j )j1k 1mmn(x j ) y ji ( x j )( x j )ckj1j1k 1mnm)i (x j )00k ( xj ) i ( xj )0(x j ) y j i ( xj )ck( xj ) k ( x j ) i ( xj ) 0j 1k 1j 1记作：n( y, i ( x)ck ( k , i )0k 1亦即（设 y( y1 , y2 ,., ym )T ）nck ( k ,i ) ( y, i ) i 1,2,., nk 1注意：上式中未知数是 ck ，是个线性方程组。腹有诗书气自华精品文档你我共享写成矩阵形式为(1 ,1 )(

31、1,2 ).(1,n )c1( y,1 )(2 ,1 )(2 ,2 ).(2 ,n )c2( y,2 )( n , 1 )( n , 2 ).( n ,n )cn( y,n )解此方程组得到最佳拟合参数，即得到表达式 Q( x)nck k ( x) 。k1三、最速下降法很多问题都可以化成求最小值的问题。例如 1：多项式拟合m 对实验数据或 yf ( x) 的一个表函数 ( xi , yi ) ， i1,2,3,., m 。n取 n 次多项式 (x)a0 a1x a2 x2 . a xnai xi，i 0其中 a0 , a1,., an 时代定系数。我们的目的就是要在 Pn ( x) 中找出一个

32、函数（确定一组 a0 , a1,., an ）n( x)aixi使偏差 yi( xi ) 的平方和i 0m2(a)(a0 , a1 ,., an )yi ( xi ) 达到极小。i 1例如 2：函数拟合用 Q( x)nck k ( x)k1去逼近函数 yf ( x) ，使二者在各点 xi 上偏差的加权平方和腹有诗书气自华精品文档你我共享m2(c)(c1 , c2 ,., cn )( xj )f (x j ) Q( xj )j 1nck k ( x) 所组成的函数在由 1 (x), 2 (x),.,n ( x) 的一切可能的线性组合k1类 H 中为极小。（确定一组 c0 , c1 ,., cn

33、）。例如 3：一般拟合m对实据 ( xi , yi , ei) ， i 1,2,3,., m用(x) 拟合数据，(x) 可以是函数，也可以是很长的一段计算过程 。( x) 中有 c0 , c1,., cn n 个参数( x)要确定一组参数 c0 , c1,., cn 使得(x) 能尽可能好的反映数据 。什么叫好呢？对计算机来说就是一个数的大小问题，所以你可以根据实际情况构造各种判断好坏的表达式。常用的有：m22m22yi( xi )加权重因子(xi)yi( xi )i 1i 1myi(xi )2m2yi( xi )2加权重因子2( x )EiiEii1i 1m 对实据 ( xi , yi ,

34、ei ) ， i1,2,3,., m 是确定的，2 的大小完全由 c0, c1,., cn 确定 。也即2是 c , c,., c的函数。2 (c , c ,., c)0 1n0 1n例如 4：方程组求解腹有诗书气自华精品文档你我共享f1 (x1, x2 ,., xn )b1fi ( x1, x2 ,., xn ) bi.f n (x1, x2 ,., xn )bnn22bifi (x1 , x2 ,., xn )*构造i 1当是解时为零，近似解不为零，越接近零近似越好。问题又归结为求一组数 x1, x2 ,., xn 使得* 尽可能小 。做法和前边的最小二乘法不一样，不是一次算出，而是逐次逼近。这类求最小值问题经常用最速下降法求。问题归结为求某函数某的最小值，不妨设函数