外观
Lesson 5 插值法 I
约 2975 字大约 10 分钟
2026-03-13
内插面对的基本问题大致如下:假定我们在一系列控制变量(自变量) x 取值处——例如我们在 x=x0,x1,⋯,xn 这 (n+1) 个两两不同的点处——获得了相应的函数 y=f(x) 的数值: y0,y1,⋯,yn,并且如果我们对函数 y=f(x) 的宏观性状做某些合理的假定(例如,它是什么类型的函数等等),我们是否可以根据这些信息完全确定,或者近似地确定这个函数本身?如果能够“近似地”确定这个函数,那么我们就不再需要去测量其他 x 处(即 x=xi,i=0,⋯,n 处)的函数值,而可以直接“近似地”计算它。这种方法是有意义的是因为一方面在我们没有测量同时又希望了解的那些点处的实验可能是十分耗时,甚至是不可能实现的;另一方面,我们有时候需要在一定的定义域内获得相应物理量 Y 的一个简单的表达式,这个表达式可能可以用于关于这个函数的进一步计算。
不失一般性,我们可以假设已经测量的这些自变量的点是按照由小到大排列的。它们被称为内插问题中的 支撑点 (support points) 或者 节点 (nodes)。即: x0<x1<⋯<xn 而整数 n 则称为支撑点的数目。我们要求内插的函数 f(x) 必须满足:
f(xi)=yi ,i=0,1,⋯,n .
显然,支撑点数目越大,我们对函数的了解就越精细。如果我们希望计算的其他点 x∈[x0,xn],这个问题称为内插问题;而如果 x<x0 或 x>xn,这个问题称为外推问题。当然,笼统地说我们也可以统称其为内插问题。
外推则不太一样。这时候我们感兴趣的是函数在某个已知区间之外的点的函数值。而且,如果我们不能够对函数的形式进行任何的限制的话,外推可能得到完全疯狂的结果。
注意内插问题与我们后面会讨论的 拟合问题 是不太一样的。内插问题中我们假定在支撑点 xi 处测得的函数值 yi 是严格的。在拟合问题中,我们一般是要考虑 yi 的误差的。但是在内插问题中,我们认为在支撑点处的值是严格的。因此,在待定函数之中的参数的数目一般也与支撑点的数目相同。所不同的是函数的形式。函数形式的选择则依赖于我们对问题的物理理解。一般来说,如果我们认为待定的函数在我们研究的区间应当是无限光滑解析的,我们会选择诸如多项式这类的函数;反之,如果我们认为所研究的函数在该区间中可能有奇异性,我们则可以选择具有极点的分式函数。在下面的几节中我们依次来讨论这些内插方法。最后我们会讨论使用样条函数进行分段拟合的情况,这种情形是函数在整个区间中并没有一个统一的形式。样条函数特别是三次样条函数在计算机图像学中有非常广泛的应用。
本章中我们还将讨论函数的计算问题。这包括两种情形。一直情形是函数具有明确的解析表达式。这个表达式可能是由级数或者一个积分表达式给出。我们需要的是从数值上准确地计算这个函数值。另外一类是首先利用某种近似来描写目标函数,然后再计算相应的函数。这一类与内插问题有类似之处。因此函数的计算无疑是我们后面几章讨论函数的积分、方程求根以及函数极值问题的基础。
Part 1 多项式内插
本节中我们讨论多项式的内插。考虑区间 [x0,xn] 上面的 n 次多项式,
Pn(x)=a0+a1x+⋯+anxn .
它具有 (n+1) 个参数: ai, i=0,1,⋯,n。我们试图利用这个多项式来解决区间 [x0,xn] 上面的内插问题。也就是说,我们希望有, Pn(xi)=yi,i=0,⋯,n。这个问题的一般解是著名的拉格朗日多项式。
我们考虑如下形式的 n 次多项式:
Lj(x)=0≤m≤n,m=j∏xj−xmx−xm ,j=0,1,⋯,n .
按照构造,这个多项式满足
Lj(xi)=δij ,0≤i,j≤n .
于是如果我们令
Pn(x)≡i=0∑nyiLi(x) ,
我们可以直接验证它满足 Pn(xi)=yi, i=0,1,⋯,n。这个公式一般称为拉格朗日内插公式 (Lagrange, 1795)。
另外一种多项式的内插方法是牛顿提出来的。它的特点是随着支撑点的增多,在利用已有的多项式的情形下,逐渐地提高多项式的阶数。
⎩⎨⎧N(x)=∑i=0naini(x) ,ni(x)=∏k=0i−1(x−xk) , n0(x)≡1 .
这个多项式的好处是,如果我们已经从 n 个支撑点的信息中获得了 ni(x),i=0,1,⋯,n−1,现在假设我们又增加了一个支撑点 (xn,yn),我们需要的仅仅是计算 an 和 nn(x)。显然, nn(x)=(x−xn)nn−1(x) 而系数 an 则可以通过求解下面的线性方程组来获得:
111⋮1⋯x1−x0x2−x0⋮(xn−x0)⋯⋯(x2−x0)(x2−x1)⋱⋯0⋮⋮⋮∏i=0n−1(xn−xi)a0⋮⋮an=y0⋮⋮yn .
这是一个下三角矩阵,因此可以利用前一章中讨论过的方法进行求解。同时,加入新的支撑点仅仅需要重新求解线性方程组的最后一行即可。
另外一种经常使用的多项式内插方法是所谓的 Neville 算法。我们下面就来介绍这种算法。假定我们的 n+1 个节点为 (xi,fi), i=0,1,⋯,n。我们利用下面的递推关系定义一系列多项式:
Pi(x)Pi0i1⋯ik(x)≡fii=0,1,2,⋯,n .=xik−xi0(x−xi0)Pi1⋯ik(x)−(x−xik)Pi0⋯ik−1(x)(a)
根据这个定义并且利用数学归纳法不难证明,这个多项式满足
Pi0i1⋯ik(xij)=fij ,j=0,1,⋯,k .
因此它满足过各个节点的条件。构造过程可以直接按照迭代式 (a) 逐层计算,例如:
例如,第一层递推中的各个函数可以根据迭代式 (a) 由相邻节点计算得到:
P01(x)P12(x)P23(x)=x1−x0(x−x0)P1(x)−(x−x1)P0(x) ,=x2−x1(x−x1)P2(x)−(x−x2)P1(x) ,=x3−x2(x−x2)P3(x)−(x−x3)P2(x) ;
然后可以进一步获得 k=2,3⋯ 的各个多项式。当我们需要增加一个节点的时候,只需要计算新增节点所对应的边界数据,已有的中间数据无需重新计算。一个常用的记号上的简化是,
Ti+k,k=Pi,i+1,⋯,i+k(x) .
采用这个记号后,递推关系可以简写为:
重要
Algorithm 4 多项式内插的 Neville 算法
Require: 给定 (n+1) 个节点信息: (xj,fj), j=0,1,⋯,n, Neville 算法给出一个 n 阶的多项式恰好过给定的 (n+1) 个节点。
【计算量】: 对于每个给定的 x 大约 O(n2).
首先填充初始值: 令 Tj0=fj, j=0,1,2⋯,n.
for k=1,⋯,j do
Tj,k=xj−xj−k(x−xj−k)Tj,k−1−(x−xj)Tj−1,k−1 ,=Tj,k−1+x−xjx−xj−k−1Tj,k−1−Tj−1,k−1 .
end for
最后的 n 阶多项式多项式由 Tn,n 给定。
对于多项式内插而言,我们有如下的误差估计。假设我们需要描写的函数 f(x) 在区间 [x0,xn] 上至少具有 (n+1) 阶的导数,那么对于任意的 x∈[x0,xn],我们一定可以找到一个 ξ∈[x0,xn] 使得,
f(x)−Pn(x)=(n+1)!ω(x)f(n+1)(ξ) ,
其中 ω(x)=(x−x0)⋯(x−xn)。也就是说,我们必定有
∣f(x)−Pn(x)∣≤(n+1)!∣xn−x0∣nx0≤ξ≤xnmax∣f(n+1)(ξ)∣
虽然我们有这个估计,但是它并不一定能够保证随着 n 的增加,误差就必定减小。
初看起来这个构造似乎已经完美解决了内插问题,至少是解决了利用进行多项式内插问题。但是实际上它是有一定的问题的。这个问题是当我们测量的支撑点的数目增加时,多项式的次数也随之增加。我们会发现,尽管在所有的支撑点处拉格朗日内插公式都满足 Pn(xi)=yi,但是在那些不是支撑点的地方,拉格朗日内插多项式可能与我们期望内插的函数——这个函数的形式我们并不清楚——可能相差很远。这个事实的另外一种体现是,当我们希望在已有的 (n+1) 个支撑点之外再增加哪怕一个支撑点,比如 (xn+1,yn+1),我们需要完全重新计算整个内插多项式的所有系数。也就是说在整个区间来看, ∣Pn(x)−Pn+1(x)∣ 可以很大,尽管在所有的支撑点处它们的差别等于零。这个现象被称为 Runge 现象 (Runge's phenomenon)。
Runge 当年考虑的函数具体的形式为
f(x)=1+25x21 .
这个函数一般称为 Runge 函数。对于 n=5 和 n=9 的拉格朗日多项式,它们分别在相应的支撑点处与 Runge 函数吻合;但是只要不在这些支撑点处,随着内插阶数的升高,拉格朗日多项式就会出现越来越严重的震荡,并可能严重偏离原函数。这就是所谓的 Runge 现象。事实上可以证明
n→∞lim(−1≤x≤+1max∣f(x)−Pn(x)∣)=+∞ .
也就是说,随着内插阶数的升高,内插多项式与原函数的最大偏离的绝对值会发散。
在相当多的应用中,内插的目的恰恰是希望得到一定区间内 没有测量 过的点处的函数值。如果我们的内插函数在支撑点之外的地方可能严重偏离实际的函数,那这种偏差对这一类的应用是不可原谅的。因此我们希望克服这种弊端。具体到 Runge 函数本身,使用分式拟合就可以很好地解决这个问题。或者使用所谓的三次样条函数来进行内插。由于 Runge 函数本身就是一个分式的形式,如果使用正确的分式内插可以完全准确地确定这个函数;而如果我们使用三次样条函数来进行内插,我们也可以获得相当不错的结果。
更新日志
2026/10/7 03:54
查看所有更新日志
98079-修正于76ef6-computational-physics于