拉格朗日多项式

在科学研究和工程实践中，我们经常需要根据若干个离散采样点来推测采样点之间的函数值。一种常见的处理方式就是，找到一个经过所有采样点的连续函数， 然后以函数值作为目标点的估计值。这个过程通常被称为插值。

解决插值的方法多种多样。我们可以把采样点 \((x_{i-1}, y_{i-1})\) 与 \((x_i, y_i)\) 用一条线段连起来，这样就可以得到一个分段函数。 分段函数虽然是连续的，但是它在采样点处通常不可导。多项式插值在数学上就很简洁美观。给定 \(n+1\) 个互不相同的采样点 \((x_0,y_0), \cdots, (x_n, y_n)\)， 需要找到一个不超过 \(n\) 次的多项式 \(P(x)\)，使其恰好经过每一个数据点，即 \(P(x_i) = y_i, i \in {0,1, \cdots, n}\)。

我们说多项式插值在数学上很优美，体现在两个方面。其一是不超过 \(n\) 次的多项式一定存在而且唯一，我们可以通过一个简单的规则构造出来。 其二是在闭区间 \([a,b]\) 上的任何连续函数，都可以用一个多项式来逼近。

1. 拉格朗日多项式

对于采样点 \((x_i, y_i)\) 我们可以构造如下的一个函数:

$$ \begin{equation}\label{f1} L_i(x) = \frac{(x - x_0)\cdots(x - x_{i-1})(x - x_{i+1}) \cdots (x - x_n)}{(x_i - x_0)\cdots(x_i - x_{i-1})(x_i - x_{i+1}) \cdots (x_i - x_n)} y_i \end{equation} $$

显然 \(L_i(x_i) = y_i\) 的，而 \(L_i(x_j) = 0, i \neq j\)。本质上 \(L_i(x)\) 是一个 \(n\) 次多项式。如果将所有 \(n+1\) 个 \(L_i(x)\) 累加起来， 得到的仍然是一个 \(n\) 次多项式 $$ \begin{equation}\label{f2} P(x) = \sum_{i = 0}^n L_i(x) = \sum_{i=0}^n y_i \prod_{i \neq k}\frac{x - x_i}{x_k - x_i} \end{equation} $$ 并且 \(P(x)\) 会经过所有采样点\((x_i, y_i)\)。这就是著名的Lagrange 多项式。 假设还有一个不超过 \(n\) 次的多项式 \(Q(x)\) 也能经过所有采样点。那么 \(P(x) - Q(x)\) 也是一个不超过 \(n\) 次的多项式。 我们知道一个 \(n\) 次多项式最多有 \(n\) 个互不相同的零点。但是对于 \(n + 1\) 个 \(x_i, y_i\) 都有 \(P(x_i) - Q(x_i) = 0\)，即 \(n+1\) 个零点。 这只能说明 \(P(x) - Q(x) \equiv 0\)，即 \(P(x) = Q(x)\)。

对于上述定理，如果 \(x = x_i, i \in \{0, \cdots, n\}\)，那么根据式 \((\ref{f2})\) 有 \(f(x) = P(x)\)。并且对于开区间 \((a,b)\) 上任意 \(\varepsilon(x_i)\)， 上式 \((\ref{f3})\) 左侧的第二项都为 0。显然对于 \(x_i\) 上述定理 2 成立。

如果 \(x \neq x_i, i \in \{0, \cdots, n\}\)，我们构造如下的一个辅助函数。由于 \(f \in C^{n+1}[a,b]\)，而多项式都是无穷阶可导 \(P \in C^{\infty}[a,b]\)， 所以该函数也是 \(n+1\) 阶可导，即 \(\varphi \in C^{n+1}[a,b]\)。

$$ \varphi(t) = f(t) - P(t) - \left[f(x) - P_n(x) \right] \prod_{i = 0}^n \frac{(t - x_i)}{(x - x_i)} $$

那么对于某个特定的 \(x\) 都有

$$ \varphi(x) = f(x) - P(x) - \left[f(x) - P_n(x) \right] \prod_{i = 0}^n \frac{(t - x_i)}{(x - x_i)} = f(x) - P(x) - \left[f(x) - P_n(x) \right] \cdot 1 = 0 $$

对于 \(t = x_i\)，有

$$ \varphi(x_i) = f(x_i) - P(x_i) - \left[ f(x) - P(x) \right] \prod_{j = 0}^n \frac{(x_i - x_j)}{(x - x_j)} = 0 - \left[f(x) - P_n(x) \right] \cdot 0 = 0 $$

所以 \(\varphi(x_i)\) 在区间 \([a,b]\) 上至少有 \(n + 2\) 个不同的零点 \(x, x_0, \cdots, x_n\)。 那么根据中值定理在 \(x, x_0, \cdots, x_n\) 的两点之间上一定存在一个数 \(\varepsilon_i\) 使得 \(\varphi'(\varepsilon_i) = 0\)。即 \(\varphi'\) 有 \(n - 1\) 个不同零点。 反复使用中值定理，可得 \(\varphi''\) 有 \(n\) 个不同零点，最终一定能找到一个数 \(\varepsilon\) 使得 \(\varphi^{(n+1)}(\varepsilon) = 0\)，有

$$ \begin{equation}\label{f4} \varphi^{(n+1)}(\varepsilon) = f^{(n+1)}(\varepsilon) - P^{(n+1)}(\varepsilon) - \left[f(x) - P(x)\right] \frac{d^{n+1}}{dt^{n+1} } \left[\prod_{i=0}^n \frac{t - x_i}{x - x_i} \right]_{t = \varepsilon} = 0 \end{equation} $$

由于 \(P(x)\) 最多有 n 次项，所以 \(P^{(n+1)}(\varepsilon) = 0\)。\(\prod_{i=0}^n \left(\frac{t - x_i}{x - x_i}\right)\) 则是一个 \(n+1\) 次的多项式， 最高次项系数为 \(\frac{1}{\prod_{i=0}^n(x-x_i)}\)。因此

$$ \frac{d^{n+1}}{dt^{n+1} } \left[\prod_{i=0}^n \frac{t - x_i}{x - x_i} \right] = \frac{(n+1)!}{\prod_{i=0}^n(x-x_i)} $$ 将上式代入式 \(\ref{f4}\) 有 $$ \begin{equation}\label{f5} f^{(n+1)}(\varepsilon) - 0 - [f(x) - P(x)] \frac{(n+1)!}{\prod_{i=0}^n(x-x_i)} = 0 \Longrightarrow f(x) = P(x) + \frac{f^{(n+1)}(\varepsilon)}{(n+1)!} \prod_{i=0}^n(x - x_i) \end{equation} $$

证毕。

式(\(\ref{f3}\))中最后的余项说明，拉格朗日多项式的逼近误差主要有两个方面的因素。 其一是 \(f^{(n+1)}(\varepsilon)\)，说明原函数的高阶导数越小，误差就越小。对应的原函数就越平滑。 其二是 \(\prod_{i=0}^n(x - x_i)\)，说明采样点的位置分布对误差也有很大影响。

2. 拉格朗日多项式的实现

如果我们严格按照拉格朗日多项式的定义 \((\ref{f2})\) 计算一个插值点 \(x\) 的值，需要两层循环才能完成，因为它有一个累加 \(\sum\) 运算和一个累乘 \(\prod\) 运算。所以其计算复杂度是 \(O(n^2)\) 的。 具体实现参见例程。 如果采样点数 \(n\) 比较多，或者需要插值的点很多，这种插值方式就会很慢。 我们之前曾提到过，对于一个 \(n\) 次多项式，通过嵌套乘法可以在 \(O(n)\) 的复杂度内完成求值计算。 所以将式 \((\ref{f2})\) 展开写出多项式的系数是有必要的。

展开式 \((\ref{f2})\) 最麻烦的地方在于，计算 \((x - x_0)\cdots(x - x_{i-1})(x - x_{i+1}) \cdots (x - x_n)\)。这本质上是 \(n\) 个一次多项式相乘。 假设目前我们已经完成了前 \(k\) 个一次多项式的乘法，将得到一个 \(k\) 多项式 \(P_k(x)\)，如下式:

$$ \underbrace{(x - x_0)\cdots(x - x_{k-1})}_{前k个}(x - x_{k}) \cdots (x - x_n) = \underbrace{(a_k x^k + \cdots + a_0)}_{P_k(x)}(x - x_{k}) \cdots (x - x_n) $$

        for (int k = 0; k < L_i.size(); ++k) {
            next_L_i[k + 1] += L_i[k];
            next_L_i[k]     -= L_i[k] * x_nodes[j];
        }

那么 \(P_k(x)\) 与下一个一次多项式 \(x - x_{k}\) 相乘，有下式。其右边可以通过右侧所示代码片段实现。

$$ \underbrace{(a_k x^k + \cdots + a_0)}_{P_k(x)}(x - x_{k}) = (a_k x^{k+1} + \cdots + a_0 x) - (a_k x_{k} x^{k} + \cdots + a_0 x_{k}) $$

完成上述多项式展开之外，我们还是需要在一个累加 \(\sum\) 运算和一个累乘 \(\prod\) 运算的双层中完成多项式所有系数的计算。 具体实现参见例程。 在工程上我们很少这样直接展开，一方面是因为 \(O(n^3)\) 的时间复杂度，更大的问题是计算误差会在多次乘加运算之后累积，导致数值不稳定。

3. Runge 现象

直觉上，应该是在闭区间 \([a,b]\) 上加密采样点，提高多项式的次数，那么得到的插值多项式 \(P(x)\) 就应该收敛到被插值的函数 \(f(x)\) 上。 事实上，并非如此。参考定理 2 中式\((\ref{f3})\) 描述的拉格朗日多项式逼近误差，也可以发现它并不能保证 \(n \to \infty\) 的时候，误差趋于零。

关于这一点，Carl Runge 举了一个例子:

$$ f(x) = \frac{1}{1 + x^2}, x \in [-5, 5] $$

在闭区间 \([-5, 5]\) 上作等距采样，当采样点数量不断增大时，得到的插值多项式在闭区间的中间位置拟合的很好，但是靠近区间两端的地方， 多项式就会剧烈上下震荡，并且振幅呈发散的趋势。这个现象被称为 龙格现象(Runge Phenomenon)。

如下图所示，我们通过例程， 分别构造了 \(5+1, 7+1, 9+1, 15+1, 17+1\) 个采样点的 Lagrange 多项式，分别对应图中 \(p5, p7, p9, p15, p17\) 的曲线。 可以明显看出 p15 和 p17 在两端震荡，而且 p17 的振幅显然大于 p15。

4. 完