第四章插值和曲线拟合

第四章插值和曲线拟合第一页，共三十一页，编辑于2023年，星期五第一节插值法的基本理论一、插值问题设函数y=f(x)给出了一组函数值yi=f(xi),i=0,1,…,n，或者给出了如下的一张表

x0,x1,x2,…,xny0,y1,y2,…,yn构造一个简单的函数φ(x)

作为f(x)的近似表达式，以满足

φ(xi)=yi

,i=0,1,…,n我们称这样的问题为插值问题。

其中φ(xi)=yi

称为插值原则;φ(x)称为f(x)的插值函数；f(x)称为被插值函数;x0,x1,x2,…,xn称为插值基点(或节点)。根据插值原则求其余点x的函数值φ(x)称为插值,x称为插值点；根据插值原则求f(x)近似函数φ(x)的方法称为插值法。第二页，共三十一页，编辑于2023年，星期五插值法的几何意义插值法的几何意义就是通过n+1个点:(xi,yi)(i=0,1,2,…,n)

作一条近似曲线y=φ(x)代替y=f(x)。如下图所示。xxnx2x1x0Xn-1y0(x0,y0)(x1,y1)(x2,y2)(xn-1,yn-1)(xn,yn)…y=φ(x)y=f(x)第三页，共三十一页，编辑于2023年，星期五插值函数φ(x)的类型在插值问题中，插值函数φ(x)的类型可有不同的选择,如代数多项式、三角多项式、有理函数等，但是最简单而常用的是代数多项式，这时就称为代数多项式插值。在本章，我们主要讨论代数多项式插值。

代数多项式插值的任务就是根据n+1个点

x0,x1,x2,…,xny0,y1,y2,…,yn构造一个次数不超过n的多项式

Pn(x)=a0+a1x+a2x2+…+anxn使满足插值原则

Pn(xi)=yi

,i=0,1,…,n。

Pn(x)称为f(x)的n次插值多项式。第四页，共三十一页，编辑于2023年，星期五二、插值多项式的误差

函数f(x)用n次插值多项式Pn(x)近似代替时，截断误差记为

Rn(x)=f(x)-Pn(x)称Rn(x)为n次插值多项式Pn(x)的余项。定理设函数f(x)在包含基点x0,x1,x2,

…,xn

的区间[a,b]上具有n+1阶导数,Pn(x)为满足Pn(xi)=yi的n次插值多项式,则对任一点x∈[a,b],总存在相应的点,使其中wn+1(x)=(x-x0)(x-x1)…(x-xn)第五页，共三十一页，编辑于2023年，星期五第二节拉格朗日插值为了得到n次的拉格朗日插值多项式，我们从最简单的一次、二次插值开始。一、一次插值（线性插值）已知x0x1求P1(x)y0y1

因P1(x0)=y0,P1(x1)=y1

所以P1(x)=y0+(y1-y0)/(x1-x0)*(x-x0)（线性插值多项式）上式可改写为：

P1(x)=y0L0(x)+y1L1(x)（拉格朗日线性插值多项式）

L0(x)=(x-x1)/(x0-x1),L1(x)=(x-x0)/(x1-x0)L0(x)、L1(x)特点：

L0(x)=1,x=x0L1(x)=1,x=x1

0,x=x1,0,x=x0第六页，共三十一页，编辑于2023年，星期五线性插值举例例已知1001/2=10，1211/2=11求1151/2解P1(x)=y0+(y1-y0)/(x1-x0)*(x-x0)P1(115)=10+(11-10)/(121-100)*(115-100)

或

P1(x)=y0L0(x)+y1L1(x)

P1(115)=10*(115-121)/(100-121)+11*(115-100)/(121-100)第七页，共三十一页，编辑于2023年，星期五

二、二次插值(抛物线插值)

二次插值问题：已知f(x)在三个互异点x0,x1,x2的函数值y0,y1,y2,要构造次数不超过二次的多项式P2(x)=a0+a1x+a2x2,使满足

P2(xi)=yi,i=0,1,2

设P2(x)=L0(x)y0+L1(x)y1+L2(x)y2,则当x=x0时,P2(x0)=y0L0(x)=1,L1(x)=0,L2(x)=0当x=x1时，P2(x1)=y1L0(x)=0,L1(x)=1,L2(x)=0当x=x2时，P2(x2)=y2L0(x)=0,L1(x)=0,L2(x)=1由上知L0(x)=1,x=x00,x=x1,x2令L0(x)=A0(x-x1)(x-x2)则A0=1/[(x0-x1)(x0-x2)]所以L0(x)=(x-x1)(x-x2)/[(x0-x1)(x0-x2)]同理可得L1(x)=(x-x0)(x-x2)/[(x1-x0)(x1-x2)]L2(x)=(x-x0)(x-x1)/[(x2-x0)(x2-x1)]综上可得P2(x)=y0(x-x1)(x-x2)/[(x0-x1)(x0-x2)]+y1(x-x0)(x-x2)/[(x1-x0)(x1-x2)]+y2(x-x0)(x-x1)/[(x2-x0)(x2-x1)]该式称为拉格朗日二次插值多项式。第八页，共三十一页，编辑于2023年，星期五二次插值举例例已知函数y=f(x)的观测数据如下表所示，试求其拉格朗日插值多项式，并计算f(1.5)的近似值。

4

-1

2

y

2

1

0

x解P2(x)=y0(x-x1)(x-x2)/[(x0-x1)(x0-x2)]+y1(x-x0)(x-x2)/[(x1-x0)(x1-x2)]+y2(x-x0)(x-x1)/[(x2-x0)(x2-x1)]=2*(x-1)(x-2)/[(0-1)(0-2)]+(-1)*(x-0)(x-2)/[(1-0)(1-2)]+4*(x-0)(x-1)/[(2-0)(2-1)]=4x2-7x+2f(1.5)≈

P2(1.5)=4*1.52-7*1.5+2=0.5第九页，共三十一页，编辑于2023年，星期五三、n次拉格朗日插值仿照P2(x)的构造方法，可得出

Pn(x)=L0(x)y0+L1(x)y1+…+Ln(x)yn其中

L0(x)=[(x-x1)(x-x2)…(x-xn)]/[(x0-x1)(x0-x2)…(x0-xn)]Lk(x)=[(x-x0)…(x-xk-1)(x-xk+1)…(x-xn)]/[(xk-x0)…(xk-xk-1)(xk-xk+1)…(xk-xn)](k=0,1,…,n)这就是n次拉格朗日插值多项式。也可写为

或

其中w(x)=(x-x0)…(x-xn)第十页，共三十一页，编辑于2023年，星期五

n次拉格朗日插值举例

例

已知函数表

x1.12751.15031.17351.972y0.11910.139540.159320.17903应用朗格拉日插值公式计算f(1.1300)的近似值。解P3(x)=L0(x)y0+L1(x)y1+L2(x)y2+L3(x)y3

=

……f(1.1300)≈

P3(1.1300)=0.1214

第十一页，共三十一页，编辑于2023年，星期五n次拉格朗日插值计算机实现

按n次拉格朗日插值公式实现第十二页，共三十一页，编辑于2023年，星期五分段插值

七、八次以上的高次插值在实际中很少采用。因为理论研究和实例都表明，插值基点增加并不能保证Pn(x)在非基点处逼近f(x)的精度得到提高,某些情况下甚至误差反而变大。所以总是对每个插值点x选择其附近的几个插值基点作低次内插(将x放在插值基点之间),或者采用分段低次插值(一次、二次插值)。为什么要选择x附近的几个插值基点？根据其中wn+1(x)=(x-x0)(x-x1)…(x-xn)第十三页，共三十一页，编辑于2023年，星期五第三节牛顿插值

拉格朗日插值多项式结构对称,使用方便,但公式不具备递推性,当需要增加基点时必须全部重新计算。因此，我们希望构造具有如下形式的插值多项式

Pn(x)=a0+a1(x-x0)+a2(x-x0)(x-x1)

+…

+an(x-x0)(x-x1)…(x-xn-1)这种形式的优点是便于改变基点数，每增加一个基点只需增加相应的一项即可(具有递推性)。为了确定出a0、a1、…、an,我们就需要讨论牛顿差商插值多项式。下面首先介绍差商的概念。第十四页，共三十一页，编辑于2023年，星期五一、差商及差商表

1.差商定义在区间[a,b]上，函数f(x)关于两点xi,xj的一阶差商定义为

f[xi,xj]=[f(xj)-f(xi)]/(xj-xi)f(x)关于三点xi,xj,xk的二阶差商定义为

f[xi,xj,xk]=(f[xj,xk]-f[xi,xj])/(xk-xi)f(x)关于k+1个点xi-k,xi-k+1,…,xi的k阶差商定义为

f[xi-k,xi-k+1,…,xi]=(f[xi-k+1,…,xi]–f[xi-k,…,xi-1])/(xi-xi-k)f(x)关于一个点xi的零阶差商定义为函数本身，即

f[xi]=f(xi)

不论几阶差商，差商均有对称性（任意改变基点的次序后其值不变）。即

f[x0,x1,…,xk]=f[]其中是x0,x1,…,xk的任一种排列。（证略）第十五页，共三十一页，编辑于2023年，星期五

2.差商表对于给定的基点及其函数值，我们可按表计算各阶差商，这样的表就叫差商表。如下:xi四阶差商一阶差商二阶差商三阶差商x0x4x3

x2x1f(x4)f(x3)f(x2)f(x1)f(x0)f[x3,x4]f[x2,x3]f[x1,x2]f[x0,x1]f[x1,x2,x3]f[x2,x3,x4]f[x0,x1,x2]f[x0,x1,x2,x3]f[x1,x2,x3,x4]f[x0,x1,x2,x3,x4]….....................零阶差商第十六页，共三十一页，编辑于2023年，星期五二、牛顿差商插值多项式

由差商定义和差商性质有f(x)=f(x0)+f[x0,x](x-x0)(f[x0,x]=[f(x)-f(x0)]/(x-x0))f[x0,x]=f[x0,x1]+f[x0,x1,x](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)=Pn(x)+Rn(x)Pn(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)Pn(x)由于满足Pn(xi)=f(xi)称作n次牛顿(差商)插值多项式。

Rn(x)=f(x)-Pn(x)=w(x)*f(n+1)(ζ)/(n+1)!

（w(x)=(x-x0)(x-x1)…(x-xn)）称为n次牛顿插值多项式的余项。第十七页，共三十一页，编辑于2023年，星期五牛顿差商插值多项式的两个特殊形式牛顿(差商)插值多项式为当n=1时

P1(x)=f(x0)+f[x0,x1](x-x0)即为线性插值。当n=2时

P2(x)=f(x0)+f[x0,x1](x-x0)+f[x0,x1,x2](x-x0)(x-x1)即为抛物线插值。可见,增加一个基点,只是增加了f[x0,x1,x2](x-x0)(x-x1)这一项。

注意,可以证明牛顿插值多项式与拉格朗日插值多项式是等价的,只不过形式不一样而已。所以,两者的截断误差是一样的。

Pn(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)第十八页，共三十一页，编辑于2023年，星期五牛顿差商插值多项式举例例已知函数y=f(x)的观测数据如下表所示，试用全部基点构造牛顿差商插值多项式，并用二次插值求f(3)的近似值。

解差商表为P4(x)=1+2(x-0)+0(x-0)(x-2)+(-1)(x-0)(x-2)(x-4)+(x-0)(x-2)(x-4)(x-5)=x4-12x3+44x2-46x+1

13

-4

9

5

1

f(x)

6

5

4

2

0

x

xi零阶一阶二阶三阶四阶

0

1

2

5

2

4

9

2

0

5

-4

-13

-5

-1

6

13

17

15

5

1第十九页，共三十一页，编辑于2023年，星期五牛顿差商插值多项式举例（续）用二次插值求f(3)时,取x0=2,x1=4,x2=5f(3)≈P2(3)=f(2)+f[2,4](3-2)+f[2,4,5](3-2)(3-4)=5+2(3-2)-5(3-2)(3-4)=5+2+5=12

xi零阶一阶二阶三阶四阶

0

1

2

5

2

4

9

2

0

5

-4

-13

-5

-1

6

13

17

15

5

1差商表为P2(x)=f(x0)+f[x0,x1](x-x0)+f[x0,x1,x2](x-x0)(x-x1)第二十页，共三十一页，编辑于2023年，星期五

牛顿差商插值多项式的计算机实现

按牛顿差商插值多项式公式实现。分两大步：

1.求各阶差商

2.按秦九韶方法求多项式的值。第二十一页，共三十一页，编辑于2023年，星期五第四节曲线拟合插值法和曲线拟合都是用来求列表函数f(x)(只知道一些点的函数称为列表函数)的近似函数φ(x)。插值法求出的近似曲线y=φ(x)要完全通过所有n+1个已知点(即要满足插值原则);而曲线拟合求出的近似曲线

y=φ(x)不要求完全通过所有n+1个已知点,只要求求得的近似曲线y=φ(x)能反映数据的基本趋势即可。曲线拟合求得的近似曲线y=φ(x)比插值法求得的近似曲线y=φ(x)更能反映客观实际。因为列表函数中的点往往是通过实验、测量或计算得来的,而实验、测量或计算得来的数据经常带有误差,如果要求所得出的曲线y=φ(x)通过所有n+1个已知点(xi,yi),就会使曲线y=φ(x)

保留着这些误差，而这是我们所不希望的。第二十二页，共三十一页，编辑于2023年，星期五一、曲线拟合问题设函数y=f(x)在n+1个互异点的观测数据为

x0,x1,…,xny0,y1,…,yn

构造函数φ(x)在包含全部基点的区间上“最好”地逼近(或靠近)f(x),这就是曲线拟合问题。如下图所示,就是使曲线y=φ(x)尽量靠近已知点(xi,yi)(i=0,1,2,…,n)。

xy0y=φ(x)×××××(x0,y0)(x1,y1)(xn,yn)(xn-1,yn-1)(x2,y2)εnεn-1ε2ε1ε0第二十三页，共三十一页，编辑于2023年，星期五

二、最小二乘法

假设y=φ(x)其中(φ(x)=a0+a1x+a2x2+…+amxm）

为给定的一组数据(xi,yi)(i=0,1,2,…,n)的拟合曲线，

则将这n+1个点代入φ(x)得以下式子

a0+a1x0+a2x02+…+amx0m≈y0a0+a1x1+a2x12+…+amx1m≈y1

……a0+a1xn+a2xn2+…+amxnm≈yn第二十四页，共三十一页，编辑于2023年，星期五最小二乘法(续)若将a0+a1x0+a2x02+…+amx0m≈y0a0+a1x1+a2x12+…+amx1m≈y1

……a0+a1xn+a2xn2+…+amxnm=yn中的“≈”换为“=”该式变为方程组a0+a1x0+a2x02+…+amx0m=y0

……a0+a1xn+a2xn2+…+amxnm≈yna0+a1x1+a2x12+…+amx1m=y1

方程组中有n+1个方程、m+1个未知数。若n+1=m+1,方程组有唯一解;若n+1<m+1,方程组有无穷多组解;若n+1>m+1,方程组无解(无精确解),称为超定方程组。第二十五页，共三十一页，编辑于2023年，星期五最小二乘法(续)

超定方程组无精确解但可求其近似解。那么解近似到什么程度才算近似解？这有不同的标准。若所求得的近似解使得误差平方和()达到最小，我们称这组近似解为最优近似解。根据误差平方和达到最小这一标准求最优近似解的方法就称为最小二乘法。下面具体解释一下什么是超定方程组的最小二乘法。

a0+a1x1+a2x12+…+amx1m=y1

……a0+a1xn+a2xn2+…+amxnm=yna0+a1x0+a2x02+…+amx0m=y0第二十六页，共三十一页，编辑于2023年，星期五超定方程组的最小二乘法

设超定方程组

根据高数知识，达到最小必须满足条件

a0+a1x1+a2x12+…+amx1m=y1a0+a1x0+a2x02+…+amx0m=y0

……a0+a1xn+a2xn2+…+amxnm=yn

(i=0,1,2,…,m)，