R语言中poly.calc与自定义实现的拉格朗日插值结果差异排查

问题场景 #

尝试为以下点集拟合拉格朗日插值多项式：

x0 <- c(576020169, 576020229, 576020296, 576020366)
y0 <- c(4816.391, 4908.896, 5011.254, 5116.875)

采用两种方法实现：

1. 使用polynom包的poly.calc函数：

library(polynom)
lagrangeInterpolation1 <- as.function(poly.calc(x0, y0))

2. 自定义拉格朗日插值函数：

lagrange <- function(x0, y0) {
    f <- function (x) {
        sum(y0 * sapply(seq_along(x0), function(j) {
            prod(x - x0[-j])/prod(x0[j] - x0[-j])
        }))
    }
    return(Vectorize(f, "x"))
}
lagrangeInterpolation2 <- lagrange(x0, y0)

两种方法结果差异极大，例如插值区间内的x=576020225：

lagrangeInterpolation1(576020225) # 2.486671e+18
lagrangeInterpolation2(576020225) # 4902.472 , 合理且正确

原因解析 #

• poly.calc的数值稳定性问题：
  polynom包的poly.calc函数是通过构造**单项式基（即$a_0 + a_1x + a_2x^2 + a_3x^{3$）**来求解插值多项式。由于你的`x0`数值量级达到$10}9$，计算高次幂时（比如$x^{3$会达到$10}{27}$级别），会超出浮点数的有效精度范围，导致多项式系数计算出现严重的浮点误差，最终代入计算后结果完全失真。

• 自定义拉格朗日函数的数值稳定性优势：
  自定义函数直接采用拉格朗日插值的基函数形式计算：每个基函数为$\prod_{k≠j} \frac{x-x_k}{x_j-x_k}$。这种形式下，分子分母都是相近大数的差值，会自动抵消掉大的基数，计算的是相对差值的乘积，避免了大数值幂运算带来的精度损失，因此数值稳定性好，结果准确。

内容的提问来源于stack exchange，提问作者Rafa