牛顿差商多项式

面对传统拉格朗日多项式的计算效率低、数值不稳定的问题，工程上经常使用牛顿差商多项式，有些资料里面也叫牛顿插值公式，来解决。和重心拉格朗日一样，牛顿差商也需要 $O(n^2)$ 的初始化，之后就可以在 $O(n)$ 的时间里完成插值运算了。只是牛顿差商的初始化过程是在构造差商表，求值的时候是通过嵌套乘法完成的。

我们已经证明经过 $n + 1$ 个互不相同的采样点的多项式 $P(x)$ 次数小于等于 $n$ 时是唯一的。所以牛顿差商多项式本质上与拉格朗日多项式是等价的，只是形式不同而已。

1. 牛顿差商多项式

给定 $n + 1$ 个互不相同的采样点 $\{x_0, \cdots, x_n\}$，及其采样值 $\{y_0, \cdots, y_n\}$，牛顿差商多项式是一个经过这些点的，不超过 $n$ 次的多项式，具有如下的形式:

$$ \begin{equation}\label{f1} P_n(x) = a_0 + a_1(x - x_0) + a_2(x - x_0)(x - x_1) + \cdots + a_n(x-x_0)\cdots(x-x_{n-1}) \end{equation} $$

如果这些点$(x_i, y_i)$是对函数 $f(x)$ 的采样，那么上式中 $a_n$ 就是 $f(x)$ 的 $n$ 阶差商(divided difference)，记作 $a_n \equiv f[x_0, x_1, \cdots, x_n]$。可以通过$n-1$阶的差商递推出来:

$$ \begin{equation}\label{f2} f[x_0, x_1, \cdots, x_n] = \frac{f[x_1, x_2, \cdots, x_n] - f[x_1, x_2, \cdots, x_{n-1}]}{x_n - x_0} \end{equation} $$

$f(x)$ 关于 $x_i$ 的零阶差商被定义为:

$$ \begin{equation}\label{f3} f[x_i] = f(x_i), i \in \{0, 1, \cdots, n\} \end{equation} $$

那么给定 $n + 1$ 个互不相同的采样点 $\{x_0, \cdots, x_n\}$，及其采样值 $\{y_0, \cdots, y_n\}$，我们可以写出如下的差商表:

$$ \begin{aligned} x_0: \quad & f[x_0] & f[x_0, x_1] = \frac{f[x_1]-f[x_0]}{x_1 - x_0} \quad & f[x_0, x_1, x_2] = \frac{f[x_1,x_2]-f[x_0,x_1]}{x_2 - x_0} & f[x_0,x_1,x_2,x_3] \quad & f[x_0,x_1,x_2,x_3,x_4] & f[x_0,x_1,x_2,x_3,x_4,x_5] \\ x_1: \quad & f[x_1] & f[x_1, x_2] = \frac{f[x_2]-f[x_1]}{x_2 - x_1} \quad & f[x_1, x_2, x_3] = \frac{f[x_2,x_3]-f[x_1,x_2]}{x_3 - x_1} & f[x_1,x_2,x_3,x_4] \quad & f[x_1,x_2,x_3,x_3,x_5] \\ x_2: \quad & f[x_2] & f[x_2, x_3] = \frac{f[x_3]-f[x_2]}{x_3 - x_2} \quad & f[x_2, x_3, x_4] = \frac{f[x_3,x_4]-f[x_2,x_3]}{x_4 - x_2} & f[x_2,x_3,x_4,x_5] \quad & \\ x_3: \quad & f[x_3] & f[x_3, x_4] = \frac{f[x_4]-f[x_3]}{x_4 - x_3} \quad & f[x_3, x_4, x_5] = \frac{f[x_4,x_5]-f[x_3,x_4]}{x_5 - x_3} \\ x_4: \quad & f[x_4] & f[x_4, x_5] = \frac{f[x_5]-f[x_4]}{x_5 - x_4} \quad \\ x_5: \quad & f[x_5] \end{aligned} $$

那么式$(\ref{f1})$ 就可以写成:

$$ \begin{equation}\label{f4} P_n(x) = f[x_0] + f[x_0,x_1](x - x_0) + \cdots + f[x_0,x_1,\cdots,x_n](x-x_0)\cdots(x-x_{n-1}) \end{equation} $$

我们定义了类 NewtonDividedDifference，其中有两个数组 mXs, mTable，分别用于保存采样点 $x_i$ 和差商表。差商表的计算如下面左侧代码片段所示，整个构造过程是在一个双层循环中完成的，时间复杂度是 $O(n^2)$。函数 InitTable 需要输入一个数组 y 与采样点 mXs 一一对应记录着各点的采样值。由于零阶差商值实际就是采样值，所以牛顿差商法并不能像重心拉格朗日插值那样一次初始化就可以应用于任何函数。

        void InitTable(std::vector<Scalar> const & y) {
            int n = mXs.size();
            mTable.assign(n, std::vector<Scalar>(n, 0.0));
            for (int i = 0; i < n; ++i)
                mTable[i][0] = y[i];
            for (int j = 1; j < n; ++j) {
                for (int i = 0; i < n - j; ++i) {
                    mTable[i][j] = (mTable[i+1][j-1] - mTable[i][j-1])
                                 / (mXs[i+j] - mXs[i]);
        }   }   }

        Scalar Evaluate(Scalar x) const
        {
            int n = mXs.size();
            Scalar result = mTable[0][n - 1]; 
            
            for (int i = n - 2; i >= 0; --i)
                result = result * (x - mXs[i]) + mTable[0][i];

            return result;
        }

完成差商表的构造之后，给定一个点 $x$，我们可以直接参考式$(\ref{f1})$以差商表的第一行为各项参数来计算插值。具体实现如上面右侧的代码片段所示。按照嵌套乘法，我们可以在 $O(n)$ 的时间复杂度内完成计算。

        void AddPoint(Scalar x, Scalar y) {
            int n = mXs.size();
            mXs.push_back(x);

            mTable.resize(n + 1);
            for (int i = 0; i <= n; ++i)
                mTable[i].resize(n + 1, 0.0);
            mTable[n][0] = y;
            for (int k = 1; k <= n; ++k) {
                mTable[n-k][k] = (mTable[n-k+1][k-1] - mTable[n-k][k-1]) 
                               / (mXs[n] - mXs[n-k]);
        }   }

2. 新增采样点

根据式 $(\ref{f4})$ 知，对于 $n + 1$ 个互不相同的采样点 $\{x_0, \cdots, x_n\}$，及其采样值 $\{y_0, \cdots, y_n\}$，经过这些点的牛顿差商多项式 $P_n(x)$ 中前 $i$ 项只与 $\{(x_k, y_k)| k ≤ i\}$ 有关。那么新增一个采样点，我们可以将式$(\ref{f4})$写成如下递推关系式:

$$ P_{n+1}(x) = P_n(x) + f[x_0,x_1,\cdots,x_n,x_{n+1}]\prod_{i=0}^n (x-x_i) $$

上式中，只需算出 $f[x_0,x_1,\cdots,x_n,x_{n+1}]$ 就可以了，这体现到差商表上，就是在下面新增一行 $f[xk,\cdots,x_{n+1}], 0 ≤ i ≤ n$ 的计算。

如右侧的代码片段所示，对应于计算 mTable[n-k][k]，它仅依赖于上一阶的差商 mTable[n-k+1][k-1] 和 mTable[n-k][k-1]，无需重新计算整个表格。更新的时间复杂度就是 $O(n)$ 的。

牛顿差商多项式

1. 牛顿差商多项式

2. 新增采样点

3. 完