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

所讨论的问题给复杂的函数

f(x)找一简单的函数p(x)如多项式、三角函数

f(x)。

等，并让其满足一定的条件，让其近似的取代原函数

或 有一数据表格，我们需要找一函数取近似的表征该表数据。

§1 拉格朗日（Lagrange）插值

在函数类中多项式具有最简单的性质。

p(x)a0a1x1a2x2a3x3...anxn

设

yf(x)在区间[a，b]连续的实函数已知在该区间上n+1个不同点xi的

f(xi)i1,2,...,n

函数值yi或 有数据表有

n1对数据 xiyii1,2,...,n 插值节点

我们需要找一个n次多项式

p(x)a0a1x1a2x2a3x3...anxn

使得在这些点上函数值等于插值节点的值。

yip(xi)

1、线性插值 已知两个点的函数值：x0y0x1y1

做一线性函数使得在两个节点上函数值为节点值。

y0p(x0)函数为：

y1p(x1)

p(x)l0(x)y0l1(x)y1xx0xx1y0y1 x0x1x1x01xjxi基函数li(x)为一次函数，且在节点上 li(xj)

0xxji几何意义：过两点做直线。按

2、抛物线插值 已知三个点的函数值：x0x变化量平均。

y0x1y1x2y2

做二次函数使得在三个节点上函数值为节点值。

y0p(x0)函数为：

y1p(x1)y2p(x2)

p(x)l0(x)y0l1(x)y1l2(x)y2xx0xx2xx0xx1xx1xx2y0y1y2 x0x1x0x2x1x0x1x2x2x0x2x11xjxi基函数li(x)为二次函数，且在节点上 li(xj)

0xxji3、拉格朗日插值 已知n+1个点的函数值：x0y0x1y1,....,xnyn ynp(xn)

做n次函数使得在n+1个节点上函数值为节点值。

y0p(x0)函数为：

y1p(x1),...,p(x)l0(x)y0l1(x)y1l2(x)y2,...,ln(x)yn

1xjxi基函数li(x)为n次函数，且在节点上 li(xj)

0xxjixxnxx1xx2...y0x0x1x0x2x0xnxx0xx2xxn...y1x1x0x1x2x1xn.....xx0xx1xxi1xxi1xxn......yixix0xix1xixi1xixi1xixn.....xx0xx2xxn1...ynxnx0xnx2xnxn1为了程序设计我们作如下推导：

j1nxxjx0xjxxjxixjny0j0j1nxxjx1xjn1y1.....ynj0jinyi.....j0xxjxnxj

i0nxxj0jiixxjjyi

程序结构是什么样子？

DIMENSION X(N),Y(N)

读入节点数据 X=.....

T1=0.0

DO 100 I=0,N T2=1.0 DO 200 J=0,N

IF(J.EQ.I) GOTO 200

xxj T2=T2* xix

j200 CONTINUE T1=T1+T2*yi

100 CONTINUE WRITE(*,*)X,T1 STOP END

拉格朗日插值的唯一性 已知n+1个点的函数值：x0y0x1y1,....,做n次函数

p(x)a123n0a1xa2xa3x...anx

使得在n+1个节点上函数值为节点值。

y0p(x0)y1p(x1),...,ynp(xn)

即：

xnyn

123na0a1x0a2x0a3x0...anx0y0123naaxaxax...ax0112131n1y1123naaxaxax...ax0122232n2y2....... 23na0a1x1axax...axn2n3nnnyn求得系数011x011x1a,a1,a2,a3,...,an则多项式确定

2x03nx0...x0线性方程组，系数矩阵对应的行列式：

x122x2x13...x1n3nx2...x201x12.......1x1n

范德蒙（Vandermonde）行列式

2xn3nxn...xn所以方程有唯一解。所以多项式唯一。

§2、牛顿插值公式 差分、差商

定义 4.1 称

f(x1)f(x0)f(x0,x1)x1x0为函数的

f(x)的一阶差商，

f(x1,x2)f(x0,x1)f(x0,x1,x2)为f(x)的二阶差商。

x2x0一般地，

f(x1,...,xn)f(x0,...,xn1)f(x0,x1,...,xn)

xnx0为

f(x)的n阶差商。为统一起见，补充定义f(x0)为零阶差商。

由差商定义，显然有

f(x0)f(x1)f(x0,x1)x0x1x1x0

f(x1,x2)f(x0,x1)f(x0,x1,x2)x2x0f(x0)f(x1)f(x2)

(x0x1)(x0x2)(x1x0)(x1x2)(x2x0)(x2x1)

如此类推，可以证明

f(x0,x1,...,xn)j0nf(xj)(xjx0)...(xjxj1)(xjxj1)...(xjxn)由此看出，差商的值与节点x0,x1,...,xn的排列次序无关，即

f(x0,x1,...,xn)f(x1,x0,...,xn)...f(x1,x2,...,xn,x0)这

种性质称为差商对称性。

定义4.2 设函数上的函数值为

yf(x)在等距节点xix0ih(i0,1,,n)f(xi)fi，其中h为常数，称做步长。称

fifi1fi 向前差分 fififi1 向后差分

fif(xih/2)f(xih/2)f分别为差分。

1i2fi1 中心差分 2f(x)在xi处以h为步长的一阶向前差分，一阶向后差分和一阶中心

符号、、分别称为向前差分算子，向后差分算子和中心差分算子。

由一阶差分的定义出发，可定义二阶差分

2fifi1fifi22fi1fi， 2fififi1fi2fi1fi2，

fif21i2f1i2fi12fifi1。

一般地刻定义n阶差分为

nfin1fi1n1fi nfin1fin1fi1

nfin1f如果令Ifii12n1f1i 2fi，Efifi1，则称I为不变算子，E为移位算子。由于

fifi1fiEfiIfi(EI)fi，

于是

EI同理可得

。

IE1，

E1/2E1/2

牛顿插值公式

有了差商的概念，前面介绍的线性插值与抛物插值可表示为

N1(x)f(xk)f(xk,xk1)(xxk)，

抛物插值

N2(x)f(xk1)f(xk1,xk)(xxk1)f(xk1,xk,xk1)(xxk1)(xxk)

事实上，按差商定义

f(x)f(x0)f(x0,x)(xx0)，

f(x0,x)f(x0,x1)f(x0,x1,x2)(xx1)，

f(x0,x1,x)f(x0,x1,x2)f(x0,x1,x2,x)(xx2)f(x0,x1,...,xn1,x)f(x0,x1,...,xn)f(x0,x1,...,xn,x)(xxn)反复用后一个式子代入前面的式子，可得

f(x)f(x0)f(x0,x1)(xx0)f(x0,x1,x2)(xx0)(xx1)f(x0,x1,...,xn)(xx0)(xx1)...(xxn1) f(x0,x1,...,xn,x)(xx0)(xx1)...(xxn)记

Nn(x)f(x0)f(x0,x1)(xx0)...f(x0,x1,...,xn)(xx0)(xx1)...(xxn1) （*）

Rn(x)f(x0,x1,...,xn,x)(xx0)(xx1)...(xxn)（**）

于是

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)!。

式（*）称作牛顿基本插值公式，它的各项系数就是函数的各阶差商，每增加一个插值节点，只需在原来的基础上多计算一项，这一性质称作牛顿插值公式的承继性。

利用牛顿插值公式时，一般分两步： 第一步：造差商表 节

函数一阶差商

二阶差商

F(x0,x1,x2) F(x1,x2,x3) F(x2,x3,x4) F(x3,x4,x5) …

三阶差商

F(x0,x1,x2,x3) F(x1,x2,x3,x4) F(x2,x3,x4,x5) …

。。。。

点 值 X0 F(x0)

X1 F(x1) F(x0,x1) X2 F(x2) F(x1,x2) X3 F(x3) F(x2,x3) X4 F(x4) F(x3,x4) X5 F(x5) F(x4,x5) … …

…

Xn F(xn) F(xn-1,xn) F(xn-2,xn-1,xn) F(xn-3,xn-2,xn-1,xn) 第二步：利用

Nn(x)f(x0)f(x0,x1)(xx0)...f(x0,x1,...,xn)(xx0)(xx1)...(xxn1)

例如：节点表 i Xi F(xi) 差商表 Xi

F（xi)

一阶

二阶

三阶

四阶

五阶

0 0.0 0.0

1 1.0 2.0

2 2.0 12.0

3 3.0 42.0

4 4.0 116.0

5 5.0 282.0

1 2 3 4

0 2 12 42

2 10 30

4 10

2

5 6

116 282

74 166

22 46

4 8

0.5 1

0.1

N5(x)02(xx0)4(xx0)(xx1)2(xx0)(xx1)(xx2)0.5(xx0)(xx1)(xx2)(xx3)0.1(xx0)(xx1)(xx2)(xx3)(xx4)即：

N5(x)02(x1)4(x1)(x2)2(x1)(x2)(x3)0.5(x1)(x2)(x3)(x4)

0.1(x1)(x2)(x3)(x4)(x5)

§3、厄密（Hermite）插值简介

高次插值逼近效果往往不理想，原因是高次函数变化复杂。 所以节点的增加并不一定能提高节点之间多项式程度。 例如： 函数

pn(x)与函数f(x)的逼近

1f(x)1x2在区间[-5，5]进行插值，我们取n1个等分节点进行

插值。会得到在节点上和函数值相等，节点之间误差很大。这种现象叫龙格（Runge）现象。

而利用分段低次多项式去逼近函数

f(x)的效果要优于光滑的高次多项式。

所以可用分段线性插值。线性插值计算简单，随节点密度增加越逼近函数，但缺点是光滑度不好（导数不连续）

分段三次厄密插值，就是要求函数值在节点上相等，函数的一阶导数也相等。 设函数为：

f(x)在n1个互异的节点xii0,1,2,...,n上函数值和导数

y0,y1,y2,...,yn进行分段插值

y'0,y'1,y'2,...,y'n

在区间[xk,xk1]上作插值，使得在两端点函数值和导数值为：

yk,yk1y'k,y'k1

四个方程，我们构造三次多项式：

23P(x)aaxaxax 30123使得P3(xk)ykP3(xk1)yk1 P3'(xk1)y'k1

P3'(xk)y'k令

P3(x)ykak(x)yk1ak1(x)y'kk(x)y'k1k1(x)

则ak(x),ak1(x),k(x),k1(x)为三次多项式且满足：

ak(xk)1,ak1(xk)0,k(xk)0,k1(xk)0

ak(xk1)0,ak1(xk1)1,k(xk1)0,k1(xk1)0 a'k(xk)0,a'k1(xk)0,'k(xk)1,'k1(xk)0 a'k(xk1)0,a'k1(xk1)0,'k(xk1)0,'k1(xk1)1

由此可求得最后ak(x),ak1(x),k(x),k1(x)

xk1xxkxkxxk1x[xk1,xk1]k0略去kn略去xxkxxk12(xx)(12xx)k1k1kkxxk12xxkak(x)()(12)xk1xkxkxk10

xxk12(xx)(xxk)xk1xxkk0略去k1kxxk12k(x)()(xxk)xkxxk1kn略去xkxk1 0x[xk1,xk1]

§4、数据拟合

问题的提出：有一组数据，数据很多，且往往有误差。需要有一函数反映这组数据。因为不同数据段误差不一样（可信度不同），同时数据多时函数中的

参数有限，所以，并不要求函数通过每一个数据点。但函数尽可能科学的反映数据。往往要确定函数中的参数，使得某种统计偏差尽可能的小。

或者说在某函数类中找一函数(x)，使得某判断量最小或最大。 往往取偏差的平方和：或[2i1m2[yi(xi)]2 最小。

i1myi(xi)2]最小，Ei为yi的误差。 Ei当函数不一样时，2,2的数值不一样，这种以函数为自变量，以数为函数值的函数关系数学上叫做泛函。

yf(x)yx2[yi(xi)]2i1m(x)2

一、用多项式进行最小二乘数据拟合

现在考察m对实验数据或yf(x)的一个表函数(xi,yi)，i1,2,3,,m。 函数类最简单的是多项式，于是我们首先取次数不超过n的代数多项式Pn(x)为H，即取HPn(x)。于是当y(x)Pn(x)时有

(x)a0a1xa2x2axnaixii0n，

其中a0,a1,an取遍所有可能的实数就得到Pn(x)的全体。

i我们的目的就是要在Pn(x)中找出一个函数(x)aix

i0n使偏差yi(xi)的平方和

(a)(a0,a1,an)yi(xi)i1m2达到极小。

i1,2,3,,m来说， 对于给定的一组数据(xi,yi)，(ai)以a0,a1,an为参数，

也就是对应Pi(xi)中一不同的(x)，(a)有不同的值，于是问题归结为选取

a0,a1,an使(a)取极小。

n次多项式由n1个 另外我们还要求nm。因为当nm时即n1m，系数所决定，而数据对只有m对，即使作出(xi)yi(i1,2,,m)才m个条件，而mn1，故(x)至少还有一个自由度，它是不确定的，但此时的(x)都使

(a)为0，当然就是极小了。所以对于nm的问题本身就没有意义。

,a的确定及最小的(a) 下面我们分几步来叙述使(a)达到极小的a01,...,an,a。 值(a01,...,an)第一步：a0,a1,...,an所满足的方程。

j(x)ax(a)(a0,a1,an)yi(xi) ji1m2nj0

,a根据微积分学，极小的a01,...,an应满足偏导数为零：

mnjxj)xk0,k0,1,2,...,n2(yiaii。 N+1个方程 ai1j0k即： 若令

mmyixiTk,ki1kxiSki1。

于是有方程组

Sajj0njkTk k0,1,2,...,n。

Sajj0njkTkSajj0njTk0Sa1Sa2S,...anST0a012n

a0S1a1S2a2S3,...anSn1T10Sa1Sa2S,...anST2a234n2,a,...,a应满足的方程。 这就是使取极小的a01n解线性方程组得到ak,k0,1,2,...,n

方程的系数矩阵为

S0S1....SnSS....S2n1Gn11 SnSn1....S2n 我们来证明系数据正对应的行列式的值不为零。

即detGn10，因为若detGn10，则齐次方程组Gn1a0,(a(a0,a1,an)T)，就有非零解a(a0,a1,an)T，从而将齐次方程组第k个方程乘以ak并对所有k从0到n作和有

0akajSkjakajxikjk0mj0k0j0i1nnnnm(akaj)xikji1mk0nj0knn(akxiajxij)i1mk0nj0n(akxjk)2k0k0

即akxjk0,k0nj1,2,...,m

意味着n次多项式，有不同的m个零点，且m>n,所以ak,k0,1,2,...,n全为零。

二、 用函数进行最小二乘曲线拟合

现在我们将前面用多项式进行拟合的情况推广到一般函数进行拟合。 设1(x),2(x),...,n(x)是n个在a,b上定义的连续函数。 所谓函数组1(x),2(x),...,n(x)是线性相关的，就是说存在n个不全为零

a,b上成立。

的实数a0,a1,...,an使a11(x)a22(x)...ann(x)0在若上述等式在a,b上除非a0,a1,...,an全为零时成立，就说1(x),2(x),...,n(x)是线性无关的。 现设1(x),2(x),...,n(x)在a,b上线性无关的。

yf(x)在xi之值为yif(xi)。

设ax1x2....xmb,T((x),(x),...,(x))引入向量记号：显然为m维列向量。 12m 所谓用函数1(x),2(x),...,n(x)去进行数据(xi,yi),i1,2,...,m的最小二乘曲线拟合，就是用线性拟合

nQ(x)ckk(x)k1

去逼近函数yf(x)，使二者在各点xi上偏差的加权平方和

(c)(c1,c2,...,cn)(xj)f(xj)Q(xj) **

j1m2在由1(x),2(x),...,n(x)的一切可能的线性组合数类H中为极小。 若在 ** 中取(x)1,式数据拟合了。

与前面要求nm一样，这里应要求n1m或nckk1nk(x)所组成的函

k(x)xk1。于是使极小的问题就变成前面多项

m，否则若nm，则

n个m维向量1,2,...,n，必然线性相关。

于是求 ** 的极小的问题就完全与求(a)的极小一样了。只不过这里是函数

1,2,...,n，比前面复杂，所以引入了新的内积写法。

为书写简化，引进“内积” (k,l)而称 (x)(x)(x)， jkjljj1mm(k,k)1/2((xj)k2(xj))1/2kj1 为m维列向量k的范数或长度。 其中(x)为a,b上给定的正连续函数，称为权函数。

(c)(c1,c2,...,cn)(xj)f(xj)Q(xj) Q(x)ckk(x)

j1m2nk1第一步：建立使(c)取极小的c所满足的方程。与前面一样，应有：

mn2(xj)yjckk(xj)i(xj)0 cj1k1ink(x)(x)0(x)ycjjkjij

j1k1m(x)(x)0(x)y(x)(x)c

jjijjkkjijj1mj1nk1mmnk(x)y(x)cjjij(xj)k(xj)i(xj)0

j1k1j1m 记作：

nk(,)0(y,i(x))cki

k1Ty(y,y,...,y)亦即（设） 12mck1nk(k,i)(y,i)i1,2,...,n

注意：上式中未知数是ck，是个线性方程组。 写成矩阵形式为

(1,1)(1,2)...(1,n)c1(y,1)(,)(,)...(,)(y,)c2222n221

(n,1)(n,2)...(n,n)cn(y,n)解此方程组得到最佳拟合参数，即得到表达式Q(x)c(x)。

kkk1n三、最速下降法

很多问题都可以化成球最小值的问题。 例如1：多项式拟合 m对实验数据或yf(x)的一个表函数(xi,yi)，i1,2,3,...,m。 取n次多项式(x)a0a1xa2x...ax其中a0,a1,...,an时代定系数。 我们的目的就是要在

n2naxii0ni，

Pn(x)中找出一个函数（确定一组a0,a1,...,an）

(x)aixi使偏差yi(xi)的平方和

i0(a)(a0,a1,...,an)yi(xi)i1m2达到极小。

例如2：函数拟合 用Q(x)c(x)

kkk1n去逼近函数yf(x)，使二者在各点xi上偏差的加权平方和

m2j1(c)(c1,c2,...,cn)(xj)f(xj)Q(xj)

在由1(x),2(x),...,n(x)的一切可能的线性组合类H中为极小。（确定一组c0,c1,...,cn）。 例如3：一般拟合 c(x)所组成的函数

kkk1nm对实据(xi,yi,ei)，i1,2,3,...,m

用(x)拟合数据，(x)可以使函数，也可以是很长的一段计算过程。

(x)中有c0,c1,...,cnn个参数

(x)

要确定一组参数c0,c1,...,cn使得(x)能尽可能好的反映数据。 什么叫好呢？

对计算机来说就是一个数的大小问题，

所以你可以根据实际情况构造各种判断好坏的表达式。 常用的有：

m2m22yi(xi)i1 加权重因子 22(xi)yi(xi)i1

2myi(xi)yi(xi)22 加权重因子 (xi)EEi1i1iim

m对实据(xi,yi,ei)，i1,2,3,...,m是确定的，

2的大小完全由c0,c1,...,cn确定。

22也即是c0,c1,...,cn的函数。(c0,c1,...,cn)

例如4：方程组求解 f1(x1,x2,...,xn)b1fi(x1,x2,...,xn)bi.......... fn(x1,x2,...,xn)bn2bifi(x1,x2,...,xn)构造

i1n2 *

当是解时为零，近似解不为零，越接近零近似越好。 问题又归结为求一组数x1,x2,....,xn使得*尽可能小。

做法和前边的最小二乘法不一样，不是一次算出，而是逐次逼近。 这类求最小值问题经常用最速下降法求。 问题归结为求某函数的最小值，不妨设函数为先取一组近似参数x1,...,xn

f(x1,...,xn)

f(x0,x1,...,xn)不为最小，记成矢量式f(X0)

但随x0,x1,...,xn 变化函数值也变化 变化最好使函数值减小，怎么变减小最快呢？ 由数学知识我们有-+的梯度矢量为

f(x0,x1,...,xn)(dfdfdfdf,,,...,) dx1dx2dx3dxn该方向就是增大最快方向，发方向就为减小最快方向。（最速下降的含义） 上式求时也可用数值求微分（后讲）。 所以求x带入

(1)ix(0)idfaidxii1,2,3,...,n

f(x0,x1,...,xn)得到新值。

在此基础上再次寻找方向，继续逼近。