第5章 插值方法
5.1 插值问题概述假设f(x)是某个表达式很复杂,甚至根本写不出来的实函数,且已知f(x)在某个区间[a,b]上的n+1个互异的点x0,x1,…,xn处的函数值f(x0),f(x1),…,f(xn),我们希望找到一个简单的函数y=P(x),使得
P(xk)=f(xk),k=0,1,…,n.
这就是插值问题。
如果我们找到了这样的函数y=P(x),我们就可以在一定范围内利用P(x)近似表示f(x),从而解决了相应的计算问题。
1.利用函数值列表来表示插值问题对于一个插值问题来说,我们的已知条件就是n+1个互异的点处的函数值.回顾高等数学中学习过的函数的表示方法,我们可用下面表1的形式列出已知的函数值,并简称为由表1给出的插值问题。
表1:插值问题的函数值列表
k
0
1
…
n
,


…




…

2.重要术语对于n+1个基点的插值问题,我们称:
f(x) 为被插值函数;
P(x)为插值函数;
x0,x1,…,xn为插值基点或插值节点;
P(xk)=f(xk),k=0,1,…,n为插值条件;
[a,b]为插值区间。
注释:对于早期的插值问题来说,f(x)通常是已知的,比如对数函数,指数函数,三角函数等这些问题现在已经不用插值法来计算了;对于现在的许多实际问题来说,我们并不知道f(x)的具体形式,所对应的函数值可能是由测量仪器或其他物理设备中直接读出来的,f(x)只是一个概念中的函数。
3.多项式插值对于n+1个基点的插值问题,如果要求插值函数是次数不超过n的多项式,记为Pn(x),则相应的问题就是多项式插值,并且把Pn(x)称为插值多项式。
实际上,我们所考虑的插值函数通常都是多项式函数或分段多项式函数。由于次数不超过n的多项式的一般形式为
Pn((x)=a0+a1x+a2x2+…+anxn (1)
所以只要确定了n+1个系数a0,a1,a2,an,我们便确定了一个插值多项式。
4.多项式插值的一般方法对于n+1个基点的多项式插值问题,我们完全可以用上一章中的办法来求插值多项式Pn(x)的系数,a0,a1,a2,an,它们可表为下面的线性方程组的解,所以多项式插值相对说来是很简单的。
 (2)
定理1:n+1级范德蒙(Vandermonde)行列式

不等于零的充要条件是诸x0,x1,…,xn两两互不相同。
这是代数学中很著名的一个定理,我们推荐大家阅读北京大学数学力学系编《高等代数》(人民教育出版社1978年第一版)pp78-79,基本方法是数学归纳法。
定理2:如果两个次数都不超过n的多项式P(x)和Q(x)在n+1个互不相同的点处的值相同,则这两个多项式恒等。
这也是代数学中很著名的定理,在北京大学数学力学系编《高等代数》(人民教育出版社1978年第一版)pp25-26中可找到定理的证明,基本思路是不恒为零的n次多项式不可能有多于n个的零点,而P(x)-Q(x)却有,所以它要么是高于n次的多项式,要么恒等于零。
5.多项式插值的基本结论由于线性方程组(2)的系数矩阵的行列式就是范德蒙行列式,再加上我们对插值基点互异的假设,上面的定理1表明,也就是说,对于给定的插值条件,存在唯一的插值多项式,它的系数就是线性方程组(2)的唯一解。
定理2进一步表明,不管什么形式的多项式,只要次数不超过n,而且满足相同的插值条件,那么它们就是恒等的。
6.重要提示,
上面的两个定理已经表明,对于多项式插值问题来说,插值多项式由插值条件唯一决定,所以我们在口语中所讲的n+1个基点的插值问题求解指的就是求出由一组插值条件所唯一决定的多项式。但是我们允许这个多项式有不同的数学形式,以便于数值计算。
在后面的两节中,我们将对相同的插值条件给出两种不同的插值多项式,即Lagrange和Newton插值多项式。所以他们应该是恒等的,名称不同只是各自的形式有所不同。
7.与曲线拟合的差异从形式上看,n+1个基点的插值问题与曲线拟合问题基本相同:都是由n+1个条件决定出一个多项式,也考虑寻找基函数的线性组合,但它们的原理,方法,以及应用领域都不相同,基本上是两个类别的问题。
曲线拟合主要应用于误差分析以及多元回归,数学原理是投影理论,方法是最小二乘法;而插值则应用于精确计算,应该说只是一种经验的方法,理论基础非常脆弱,不过我们可以通过验算来证实结果的有效性。
5.2 拉格朗日插值公式拉格朗日插值公式的构造对于表1所确定的n+1个基点的插值问题,相应的拉格朗日插值公式,我们记为Ln(x),实际就是把插值多项式Pn(x)表为基函数{lk(x)|k=0,1,…,n}的线性组合,其中
 (1)
而且组合系数实际就是诸y0,y1,…,yn,因此我们有
 (2)
接下来我们要证明Ln(x)的确是插值多项式,也就是满足n+1个插值条件。
证明Ln(x)的确是插值多项式注意到所有的lk(x)都是n次数多项式,所以Ln(x)是次数不超过n的多项式。再注意到

所以我们有

所以Ln(x)满足插值条件,即Ln(x)的确是插值多项式。
3.基函数的性质定理:对于给定的函数值列表
,


…




…

以及由(1)定义的基函数{lk(x),k=0,1,…,n},我们有

证明:考虑被插函数为f(x) ≡1的特殊的多项式插值问题,对于给定的插值基点x0,x1,…,xn,对应的函数值y0,y1,…,yn均为1,从而相应的拉格朗日插值多项式为

它满足插值条件

而f(x)≡1也是次数不超过n的多项式,也满足插值条件
f(xk)=yk,k=0,1,…,n
由插值多项式的唯一性立即推出我们的结论。
4.计算表格与算法说明应该说,对于给定的x,求拉格朗日插值多项式的值还是很简单的,直接按公式计算就成,所以没有太大的价值来研究算法问题。如果用计算器来计算,我们建议按如下的表格计算。
k
X[k]
Y[k]
l[k]
L[k]
0
X[0]
Y[0]
l[0]
L[0]
1
X[1]
Y[1]
l[1]
L[2]
…
…
…
…
…
n
X[n]
Y[n]
l[n]
L[n]
x
1
首先依次计算基函数l[0],l[1],…,l[n],l[k]存放基函数lk(x)的值。接下来再依次计算
L[0]=Y[0]*l[0],
L[k]= Y[k]*l[k]+L[k-1],k=1,2,…,n.
5.计算基函数值的方法拉格朗日多项式插值的计算主要是计算基函数的值,作为练习编程,编写求单个的基函数值的程序还是很有意义的。下面是求基函数的程序段:
float GetLK(float x,int k)
{ int i = 0;
float y=1.0;
for(i=0;i<=n;i++)
{ if(i==k) continue;
y *= x-x[i];
y /= x[k]-x[I];
}
return y;
}
6.补充说明拉格朗日插值虽然算法很简单,但是当n较大时,用计算器来算还是有点累赘,大家一定不要再去寻找捷径,唯一的笨办法就是老老实实地一个一个地算,把结果填在表中。可能用纸和笔作减法运算,用计算器作乘法和除法计算这种组合模式最佳。基函数值计算完毕后,可作一次检验,看累加起来是否为1,如果是,则计算正确;如果不是,则前面的计算肯定有误,应当重新计算。
基函数的等价形式记 
则有 
记 
则有 

所以有 
8.插值函数的等价形式利用lk(x)的等价形式,我们可以进一步得到Ln(x)的等价形式:

结论:拉格朗日插值多项式的xn项(首项)的系数为

在有些场合下,这些公式也能发挥作用。
5.3牛顿插值公式具有承袭性的插值多项式承袭性的含义假设对于具有k+1(k≥0)个插值基点的多项式插值问题已经得到了相应的插值多项式,如果它不能满足我们的精度要求,则希望通过增加一个插值基点,并利用已得到的来构造新的,具有k+2个插值基点的插值多项式.为此,把表为与一个修正项,记为,的和的形式是合理的,即
=+ (1)
从而只要确定了,我们即可写出.此时我们称上面的(1)是具有承袭性的插值多项式。
具有承袭性的插值多项式的一般形式注意到(1)式中的和均应满足前面的k+1个插值条件
== j=0,1,…,k
所以我们可以立即得到
=0 j=0,1,…,k (2)
再注意到在一般情况下应当是k+1次多项式,因此可表为
=… (3)
其中为待定常数.事实上它可以由最后新增的一个插值条件来确定,但我们现在不在意。

=… (4)
则有
=,
而且还有
=0 j=0,1,…,k (5)
为方便计,我们特别约定≡1,从而还有:
= (6)
对于由表1给出的一般的具有n+1个插值基点的多项式插值问题,令k分别取0,1,…,n-1,上面的(1)到(6)式仍然成立。现在我们令k依次取n-1,n-2,…,0,反复利用(1)式和(3)式即可得到
 (7)


上面的(7)式是具有承袭性的插值多项式的一般形式,它可以看作是基函数的线性组合.
确定组合系数为了使(7)式真正成为插值多项式,还需要根据插值条件确定组合系数,,…,.
在(7)式中分别令取,,…,,并利用插值条件

和(5)式,我们有
 (8)
这是一个下三角形的线性方程组,直接求它的数值解也很容易.
2.差商的定义与性质为了深入研究线性方程组(8)的解的特性,我们把计算方法教材中关于差商的定义修订如下:
设是定义在某区间内取连续变量值或离散变量值的实值函数,且已知在n+1个互异的离散点,,…,处的函数值,,…,,称
= k=0,1,…,n (9)
为在处的零阶差商;
  (10)
为在,处的一阶差商;
 
为在,,处的二阶差商.
一般地,如果定义了的k-1阶差商,则可以定义的k阶差商为:
 (11)
这里,k可取值 2,…,n.
定理1,设是定义在某区间内取连续变量值或离散变量值的实值函数,对该区间内任意k个互异的离散点以及任意异于的x,我们有
 (12)
证明:当k=0时,由

立即推得

这也就是

所以当k=0时结论成立.
假设当k取某个时结论也成立,即
 (13)
则当k=j+1时,由差商的定义得

整理后得

再代入到(13)式中,并利用

即可推得
所以结论对k=j+1也成立。
所以,(12)式对所有的都成立。证毕牛顿插值公式现在,我们只要把差商的定义以及上面的定理1与表1所确定的n+1个基点的插值问题联系起来,即可导出具有承袭性的牛顿插值多项式。
定理2:对于表1所确定的n+1个基点的多项式插值问题,对任意k=0,1,…,n-1,我们有
(14)
证明:在(12)式中令立即得证。
定理3:对于表1所确定的n+1个基点的多项式插值问题,我们有

为线性方程组(8)的解。
这个定理实际上是上面的定理2的推论。但我们还是鼓励读者根据差商的定义直接用数学归纳法证明,难度不算太大。
事实上,只要按k=0,1,…,n-1的顺序把上面的(14)式重写一遍,在最前面补上,再与(8)式对比,利用解的唯一性以及即可完成证明。
定理2和定理3表明,对于
 (15)
是相应于插值条件

的、具有承袭性的插值多项式。
特别地在(12)式和上面的(15)式中取k=n,并令
 (16)
我们有

它们分别称为牛顿插值多项式、插值余项和牛顿插值公式。
至此,我们已经得到了具有承袭性的插值多项式.
4.计算表格与计算方法说明对于给定的n+1个基点的插值问题以及自变量的任意取值x,我们建议按如下的表格计算牛顿插值多项式的系数以及多项式的值。
K
Y[k]
X[k]
F[k]
P[k]
N[k]
0
Y[0]
X[0]
F[0]
P[0]
N[0]
1
Y[1]
X[1]
F[1]
P[1]
N[1]
…
…
…
…
…
…
N
Y[n]
X[n]
F[n]
P[n]
N[n]
表中X[k]列和y[k]列为函数值列表,
F[k]存放,
P[0]=1,
P[k]=(x-xk-1)*P{k-1},k=1,2,…,n
N[0]=F[0]*P[0]
N[k]=F[k]*P[k]+N[k-1],k=,2,…,n
即使利用计算器来算,按我们这里给出的定义计算差商也更方便些,计算F[k](这一列)的程序段如下:
float GetF(int k)
{ int i;
float w;
F[k]=Y[k];
for(i=1;i<k;i++) F[k]=(F[k]-F[I])/(X[k]-X[i]);
return;
}
提示:在循环过程中,F[k]存放的实际是f[x0,x1,…,xi,xk],所以当循环到i=k-1时,所得即为f[x0,x1,…,xk-1,xk]
5.4关于多项式插值余项的讨论插值余项的定义对于一般的n+1个基点的多项式插值问题,我们在前面的三节中分别给出了三种不同的方法:一般方法、拉格朗日插值法和牛顿插值法。由唯一性定理可知,不同的方法得到的不同形式的不超过n次的多项式由于在n+1个不同的点处的值相同,所以它们实际上是恒等的,所以本节关于插值余项讨论对前面三节均有效。
定义:对于一般的n+1个基点的多项式插值问题,设f(x)为被插值函数,Pn(x)为相应的插值多项式,记Rn(x)为f(x)与Pn(x)的差,即
f(x)= Pn(x)+ Rn(x) (1)
则Rn(x)就是用Pn(x)近似替代f(x)的误差,我们称它为插值余项.
显然,由插值多项式的唯一性可以导出插值余项的唯一性。
插值余项的估计定理1:设f(x)是定义在[a,b]上的具有直至n+1阶连续导数的被插值函数,又假设n+1个插值基点满足a≤x0<x1<…<xn≤b,Pn(x)为相应的插值多项式,Rn(x)为插值余项,则对任意固定的x∈(a,b),我们有
 (2)
证明提示:把x看作是(任意)固定的点,作辅助函数

则g(t)有n+2个零点,n+1次应用罗尔定理,则存在ξ使得g(n+1)(ξ)=0,注意到x是常数,Pn(t),πn+1(t)是多项式,所以推出我们的结果并不难。
差商与导数的关系在牛顿插值公式一节中,我们曾给出了多项式插值余项的数学形式。对(a,b)内任意异于x0,x1,…,xn的x,我们有
 (3)
把它与上面的(2)式对比即可看出
 (4)
4多项式的差商利用牛顿差值公式以及余项的数学形式,我们还可以推出许多有意义的结果,比如:
任何n次多项式的n+1阶差商均为零。
由n次多项式的n+1个互异的基点的牛顿插值公式的余项Rn(x)=0,立即推得f[x0,x1,…,xn,x]=0.
这个结论也可以由上面的(4)式立即推出。
任何n次多项式的n阶差商均为常数。
考察n次多项式P(x)的任意的n+1个基点的牛顿插值多项式Nn(x),显然它们恒等,从而首项系数相同,所以P[x0,x1,…,xn]为常数。
差商具有对称性,也就是说,如果z0,z1,…,zn是x0,x1,…,xn的任意一个排列,那么f[z0,z1,…,zn]= f[x0,x1,…,xn].
证明:仅改变插值基点的排列顺序所得到的两个牛顿插值多项式是恒等的,所以他们的最高次项的系数也相等,得证。
5.6 Hermit插值与分段插值
1.多项式插值的龙格现象设被插函数为f(x)=1/(1+x2),插值区间为[-5,5],把插值区间n等分,插值基点为区间的等分点。我们会发现,随着n的增大,插值多项式灾区尖端点附近发生激烈的震荡,这就是龙格现象。
克服的办法:减小插值区间,降低多项式的次数,要求插值函数与被插函数在基点处的导数值也相同。
2.Hermit插值问题:已知f(x)在x0,x1,…,xn处的函数值y0,y1,…,yn和导数值z0,z1,…,zn,我们可以构造一个2n+1次多项式H(x),满足这2n+2个条件,即

仿照拉格朗日插值公式的构造方法,我们也可以把H(x)表为基函数αk(x),k=0,1,…,n和βk(x),k=0,1,…,n的线性组合,且组合系数相应地分别为y0,y1,…,yn和z0,z1,…,zn,则有:

显然,H(x)是2n+1次多项式,它由2n+2个常数确定,而则2n+2个常数正好由2n+2个插值条件确定。
不难看出,对任意给定的j=0,1,…,n,

可以证明,对任意给定的j=0,1,…,n及k=0,1,…,n

所以H(x)的确满值插值条件。
提示:由于诸αk(x)和βk(x)都是连乘积的形式,所以用对数求导法来计算有关的导数值更方便些。