重心插值

重心插值（barycentric interpolation）用一组只依赖节点的权重，计算通过给定数据点的插值多项式。它是拉格朗日插值多项式的一种代数表示，不是把两个邻点作线性加权，也不是用点在几何单纯形中的重心坐标作任意平均。固定节点后，每次求值只需线性数量的运算。

取三个教学数据点 $(-1,1),(0,0),(2,4)$。显然多项式 $p(x)=x^2$ 通过它们；重心公式将在不先求幂系数的情况下重新算出 $p(1)=1$，并说明节点数据改变后结果如何变化。

节点权重与两个等价公式

设节点 $x_0,\dots ,x_n$ 两两不同，数据值为 $y_0,\dots ,y_n$，令 $$\omega (x)=\prod _{j=0}^n(x-x_j),w_i=\frac{1}{\prod _{j\ne i}(x_i-x_j)}=\frac{1}{{\omega }^{\prime }(x_i)}.$$ 节点不同保证权重分母非零。拉格朗日基函数可写为 $$L_i(x)=\prod _{j\ne i}\frac{x-x_j}{x_i-x_j}=\omega (x)\frac{w_i}{x-x_i}(x\ne x_i).$$ 因此插值多项式有第一重心形式 $p(x)=\omega (x)\sum _i{w}_i{y}_i/(x-x_i)$。进一步约去公共因子，得到更常用的第二重心公式： $$p(x)=\frac{\sum _{i=0}^n\frac{w_i{y}_i}{x-x_i}}{\sum _{i=0}^n\frac{w_i}{x-x_i}},x\notin \{x_0,\dots ,x_n\}.$$ 在节点处直接定义 $p(x_i)=y_i$，避免计算形如除以零的中间项。

有理外形背后的多项式证明

先看拉格朗日表达式 $\sum _i{y}_i{L}_i(x)$，它次数至多 $n$，且在节点 $x_j$ 处只有第 $j$ 项留下，值为 $y_j$。若两个次数至多 $n$ 的多项式都匹配 $n+1$ 个不同节点，它们之差具有 $n+1$ 个不同零点，只能为零，所以插值多项式唯一。

将全部数据值暂取为一，常数多项式一已满足插值条件；由唯一性得到 $\sum _i{L}_i(x)=1$。在非节点处代入上一节的基函数表达式，即得 $$\omega (x)\sum _i\frac{w_i}{x-x_i}=1.$$ 所以分母恰为 $1/\omega (x)$，数学上非零；拿第一重心形式除以这个恒等式，便得到第二形式。虽然它写成分式，在标准权重下依旧是同一个多项式。若任意另选权重，这个保证可能消失，得到的可以是真正的有理插值函数。

三点手算

主例的权重为 $$w_0=\frac{1}{(-1-0)(-1-2)}=\frac{1}{3},w_1=\frac{1}{(0+1)(0-2)}=-\frac{1}{2},w_2=\frac{1}{(2+1)(2-0)}=\frac{1}{6}.$$ 在 $x=1$ 处，分母为 $\frac{1}{6}-\frac{1}{2}-\frac{1}{6}=-\frac{1}{2}$，分子为 $\frac{1}{6}+0-\frac{4}{6}=-\frac{1}{2}$，所以商为一。权重有负号，并非概率权重；分母也不必为正。

权重由节点决定，查询点不需要重新求多项式的幂系数。

全部权重同时乘一个非零常数，分子和分母同倍变化，结果不变。例如可把本例权重乘六写为 $(2,-3,1)$。这提供了数值缩放的自由，但不允许分别任意改变各个权重。

变式。 只把第三数据值从四改为五。插值多项式增加对应的一个基函数 $$L_2(x)=\frac{(x+1)x}{(2+1)(2-0)}=\frac{x(x+1)}{6}.$$ 所以新多项式是 $x^2+x(x+1)/6$，在一处为 $4/3$。节点未变，权重完全不用重算。这个变式也说明，某个数据值误差在查询点被放大多少，由相应基函数决定。

可运行的Python求值器

函数接受相同长度的有限实数节点与数据，节点必须不同，返回权重元组和一个求值函数。输入不被改写。这里用朴素乘积计算权重，适合规模较小、尺度适中的教学数据；不对极端节点范围下的浮点溢出作额外处理。

Python 3

from math import fsum, isfinite

def barycentric_interpolator(nodes, values):
    nodes = tuple(float(t) for t in nodes)
    values = tuple(float(y) for y in values)
    n = len(nodes)
    if n == 0 or len(values) != n:
        raise ValueError("nonempty equal-length data required")
    if not all(isfinite(z) for z in nodes + values):
        raise ValueError("finite data required")
    if len(set(nodes)) != n:
        raise ValueError("nodes must be distinct")
    weights = []
    for i in range(n):
        product = 1.0
        for j in range(n):
            if i != j:
                product *= nodes[i] - nodes[j]
        weights.append(1.0 / product)
    weights = tuple(weights)
    def evaluate(x):
        x = float(x)
        if not isfinite(x):
            raise ValueError("finite query required")
        for node, value in zip(nodes, values):
            if x == node:
                return value
        terms = [w / (x - node) for w, node in zip(weights, nodes)]
        denominator = fsum(terms)
        if denominator == 0:
            raise ArithmeticError("numerical cancellation in denominator")
        return fsum(t * y for t, y in zip(terms, values)) / denominator
    return weights, evaluate

weights, p = barycentric_interpolator([-1, 0, 2], [1, 0, 4])
assert abs(p(1) - 1) < 1e-14
assert p(0) == 0
_, changed = barycentric_interpolator([-1, 0, 2], [1, 0, 5])
assert abs(changed(1) - 4 / 3) < 1e-14
_, constant = barycentric_interpolator([2], [7])
assert abs(constant(-10) - 7) < 1e-14
try:
    barycentric_interpolator([0, 0], [1, 2])
except ValueError:
    pass
else:
    raise AssertionError("repeated nodes accepted")
print(weights, p(1))

权重的双重循环需 $O(n^2)$ 时间，每次查询的节点检查与求和需 $O(n)$，保存数据和临时项需 $O(n)$ 空间；这里复杂度中的 $n$ 指数据点数。计算上的快速表示并不保证插值问题本身条件良好。高次等距插值仍可能出现龙格现象，远离节点区间的外推也可能放大误差；应把表示的计算稳定性与节点选择导致的逼近性质分开考察。