函数插值与插值多项式

函数插值在物理学,天文学中有广泛的应用。函数差值帮助我们从一群实验散点中预测未来的数值趋势和走向。我们在这里主要讨论两个最简单的多项式插值方法,即拉格朗日插值和牛顿插值。

函数插值与插值多项式

插值法

对于函数$f(x)$,已知一些列的函数值,$(x_i, y_i)$,绘制一条穿过这些函数值的曲线,$P(x)$,使得$ P(x_i) = f(x_i) = y_i $,称函数$P(x)$为函数$f(x)$的插值函数。

在现实情况中,往往原函数$f(x)$的解析形式为知,或者过于复杂,只能通过实验的方式取得离散的函数值,即$(x_i, y_i)$,通过这些离散的散点,利用函数插值,还原出原函数的近似函数。

特别的,插值函数$P(x)$,我们考虑一个简单的形式,即多项式函数,此时,称插值函数$P(x)$为原函数的多项式插值。

假设多项式插值函数$P(x)$为一个$n$次函数,那么多项式$P(x)$即有$n+1$个系数待定,即$a_0,a_1,...,a_n$。利用待定系数法可求得,至少需要$n+1$个散点,得到一个关于系数的线形方程组,即可求出$n+1$个系数。

容易得知,对于$n+1$个线性无关的系数,其解唯一,即相对应的多项式插值函数$P(x)$唯一。

多项式原函数与多项式插值函数

容易得知,如果原函数$f(x)$是一个多项式函数,且最高次数为$n$,那么,根据对应的$n+1$个散点求的的多项式插值函数$P(x)$,则一定和原函数相同。

%matplotlib inline

import numpy as np
import matplotlib.pyplot as plt
from sympy import *
from sympy.utilities.lambdify import lambdify

x = symbols('x')
y = x * x
f = lambdify(x, y, 'numpy')

plt.axis([-5, 5, -2, 20])

xvalues = np.arange(-5, 5, 0.1)
yvalues = f(xvalues)
plt.plot(xvalues, yvalues, label='f(x)')

xvalues = np.array([1, 2, 3])
yvalues = f(xvalues)
plt.scatter(xvalues, yvalues, color='red')

plt.show()

如上图所示,对于一个二次函数$f(x)$,其函数上的任意三个点的多项式插值函数一定等于原函数。

from scipy.interpolate import lagrange

y = sin(x)
f = lambdify(x, y, 'numpy')

xvalues = np.arange(-5, 5, 0.1)
yvalues = f(xvalues)
plt.plot(xvalues, yvalues, label='f(x)')

xvalues = np.array([1, 2, 3])
yvalues = f(xvalues)
plt.scatter(xvalues, yvalues, color='red')

p = lagrange(xvalues, yvalues)
xvalues = np.arange(0, 4, 0.1)
yvalues = p(xvalues)
plt.plot(xvalues, yvalues)

plt.show()

对于上图表示的函数,原函数为一个正弦函数$sin(x)$,对于在正弦函数上的散点,利用拉格朗日二次插值,橙色函数为二次插值函数$P(x)$,可以看出,在散点附近,原函数和插值函数接近,随着越来越远离采样点,插值函数将越发偏离原函数。

误差,插值余项

对于原函数和插值函数,由于插值函数是原函数在区间$[a, b]$上的一个近似,我们即可定义插值余项$R(x) = f(x) - P(x)$即为在区间$[a, b]$上的误差。

对于一个多项式插值,由于插值函数由散点序列唯一确定,误差也即由散点序列唯一确定。

由于多项式的系数由系数的矩阵唯一确定,然而求解一个$n$元的线形方程的解将消耗非常多的时间,据此,拉格朗日和牛顿分别提出了两种求解方式。注意:两种方式只是符号形式的不同,其解的值,和误差均相等。

拉格朗日插值

对于$n+1$个散点,可以直接确定一个$n$次多项式的系数。拉格朗日插值法构造了$n+1$个$n$次多项式,并确保这$n+1$个多项式的线性组合一定穿过$n+1$个散点,由此即可确定,其线形组合的函数即为插值函数。

由于对于任何一个$x_i$,需要满足$P(x_i) = f(x_i) = y_i$。考虑一个特殊的函数集合,其中:

$$f_i(x) = \begin{cases}1, x=x_i \cr 0, x=x_j, j \neq i\end{cases}$$

其中$f_i$是一个$n$次多项式。对于$x_i$这个散点,函数值为$1$,对于其他的所有散点,函数值为$0$。即该多项式穿过了$(x_0, 0), ..., (x_i, 1), ..., (x_n, 0)$。

通过如此构造,我们即可通过线性组合得到插值函数$P(x)$:

$$P(x) = \sum_{i=0}^{n}{y_i*f_i}$$

现在构造$f_i(x)$,由于只有在$x_i$处不是零点,那么显然

$$f_i(x) = (x-x_0)...(x-x_{i-1})(x-x_{i+1})...(x-x_n)a_i$$

又由于$f_i(x_i) = 1$,那么很显然可以得到

$$f_i(x) = \frac{(x-x_0)...(x-x_{i-1})(x-x_{i+1})...(x-x_n)}{(x_i-x_0)...(x_i-x_{i-1})(x_i-x_{i+1})...(x_i-x_n)}$$

我们可以定义个一个连乘函数$w_n(x) = (x-x_0)(x-x_1)...(x-x_n)$,根据求导乘法法则,有$(fg)' = f'g + fg'$,可以从$w_n(x)$提取出$x-x_k$项,并利用乘法法则求导,则有

$$w_n'(x_k) = (x_k-x_0)...(x_k-x_{k-1})(x_k-x_{k+1})...(x_k-x_n)$$

因此,函数$P(x)$可表示为

$$P(x) = \sum_{i=0}^{n}{\frac{y_i*w_n(x)}{(x-x_i)w_n'(x_i)}}$$

牛顿插值

拉格朗日插值由一组函数的线性组合得到插值函数$P(x)$,其构造简单直观,对于n较小时,能极其有效的求解插值函数。然而,拉格朗日插值中的系数需要所有的散点联合求得,因此在任何一个散点发生变化的时候,或者增加新数据的时候,拉格朗日插值中的所有系数都需要重新计算,为此,牛顿插值提供了更好的解决方式。

首先考虑两个采样点,那么插值函数即为一次函数。一次函数解析式可直接由点写斜式给出

$$P_1(x) = \frac{y_1-y_0}{x_1-x_0}(x-x_0) + y_0$$

很显然从上面的式子可知,$P_1(x_0) = y_0, P_1(x_1) = y_1$。现在我们新增一个点$(x_2, y_2)$,由于新的插值函数依然穿过了之前的两点,那么我们考虑新的函数

$$P_2(x) = P_1(x) + r_2(x)$$

很显然,当$x = x_0, x_1$时,$r_2(x) = 0$,并且$r_2(x)$也是一个多项式,因此,可以构造

$$r_2(x) = a_2(x-x_0)(x-x_1)$$

至于$a_2$的值,我们需要利用$(x_2, y_2)$求解。于是

$$y_2 = P_1(x_2) + a_2(x_2-x_0)(x_2-x_1)$$

$$a_2 = \frac{y_2 - P_1(x_2)}{(x_2-x_0)(x_2-x_1)}$$

更一般的,我们可以得到一个$P_n(x)$的递推公式

$$P_n(x) = \begin{cases}P_{n-1}(x) + a_n\prod_{i=0}^{n-1}{(x-x_i)}, n \ge 1 \cr y_0, n = 0\end{cases}$$

如果我们将递推公式的左右两侧进行连加,那么即可消去中间项,令$a_0=y_0$,于是得到

$$P_n(x) = a_0 + \sum_{i=1}^{n}{(a_i\prod_{j=0}^{i-1}{(x-x_j)})}$$

现在我们来确定系数$a_n$。首先我们将$P_n(x)$,根据递推式,拆解为

$$P_n(x) = a_nA_n(x) + P_{n-1}(x)$$

根据定义可以知道,$P_{n-1}$为$n-1$次多项式,$A_n$为$n$次多项式。我们可以将$P_n$认为一个在高维空间的直线,那么$a_n$即为该直线的斜率。因此,$a_n$可以使用求导的形式得到,又因为$P_n(x)$是一个最高次为$n$的多项式函数,根据求导定义,通过$n$次求导以后,剩下的即为最高次的系数$a_n$。因此

$$P_n^{(n)}(x) = a_n * n!$$

遗憾的是$P_n^{(n)}(x)$是一个为知量,如此构造并不能为我们提供$a_n$的解(实际上,上面式子通过一定变换,即为牛顿插值的泰勒形式)。然而,我们可以通过求导的形式,找到求解的思路。由于每一次求导实际上是求了一个差商,而通过一次差商,我们即可消去一个低次项。例如

$$y = ax^2+bx+c$$

其差商

$$\frac{a(x_1^2-x_0^2) + b(x_1-x_0)}{x_1-x_0} = a(x_1+x_0) + b$$

由于我们的目标是求出最高次系数$a$,那么我们还需要消去方程中的$b$项。由于$b$项现在已经是低次项,同样的,我们可以对该差商,构造一个新的差商。

$$\frac{a(x_2^2-x_0^2) + b(x_2-x_0)}{x_2-x_0} = a(x_2+x_0) + b$$

由此看出,我们只需要将两次相减,取差商,即可消去$b$项,同时得到$a$项

$$a = \frac{(a(x_2-x_0) + b) - (a(x_1-x_0) + b)}{x_2-x_1}$$

均差

我们看到,对于离散形式的函数,在无法求导时,我们可以用多次差商来近似地取得函数的高阶导数。于是我们引入均差的定义

给定函数$f$和一组数据点,$(x_0, y_0), ..., (x_n, y_n)$,引入均差记号$f[x_0,...,x_n]$,

定义$f[x_i] = y_i$,并且

$$f[x_i,...,x_{i+j}] = \frac{f[x_{i+1},...,x_{i+j}] - f[x_i,...,x_{i+j-1}]}{x_{i+j}-x_i}$$

很显然的,结合上例关于二次函数差商的推导,通过均差运算,最后得到的是最高次的系数,即$a$。于是牛顿插值即可表示成

$$ P_n(x) = a_0 + \sum_{i=1}^{n}{(a_i\prod_{j=0}^{i-1}{(x-x_j)})} = f[x_0] + \sum_{i=1}^{n}{(f[x_0,...,x_i]\prod_{j=0}^{i-1}{(x-x_j)})} $$

虽然均差涉及到一个递归式的求解,但由于均差仅仅依赖于前一项,由此,我们可以利用动态规划,维护一张均差表,由前项均差推导出后项均差。