拉格朗日插值算法低阶正常但59阶计算结果异常求助

问题描述 #

参考资料实现了Lagrange插值算法，代码如下：

import numpy as np
def lagr(P) :

    A = list()
    for i in range(0, len(P)) :
        A.append(P[i][0])

    B = list()
    for i in range(0, len(P)) :
        B.append(P[i][1])

    F = []
    for i in range(0, len(A)) :
        AA = list()
        for k in range(0, len(P)) :
            AA.append(A[k])

        del AA[i]

        ff = np.poly1d(AA, True)

        aa = 1
        for j in range(0, len(AA)) :
            aa  *= A[i] - AA[j]

        f = ff * B[i] / aa
        F.append(f)

    addf = 0
    for k in F :
        addf = addf + k

    return addf

测试3个点生成二次函数时结果正常：

P = [[0, 1], [1, 4], [3, 2]]
lagr(P)
>>>poly1d([-1.33333333,  4.33333333,  1.        ])

但处理60个点生成59阶多项式时结果异常，测试代码如下：

A = [27, 35.1, 26.3, 32.5, 29.1, 32.9, 28.4, 33.5, 26.4, 32.1, 26, 31.4, 25.7, 30.9, 24, 27.8, 25, 27.6, 24.8, 30.5, 23.9, 29.6,
     27.3, 31.8, 26.1, 30.4, 25.7, 29.5, 26.4, 33, 26.4, 34, 26.9, 32.1, 28, 34.9, 26.5, 34.1, 27, 32.7, 27.2, 31.1, 27.6, 33.8,
     27.9, 34.8, 28.4, 34.2, 27.6, 33.3, 28, 32.8, 28.6, 33.1, 27.9, 32.5, 26.3, 29, 25.9, 31]

B = list(range(1, 61))

x = len(A)
H = []
for i in range(0, x) :
    h = [B[i], A[i]]
    H.append(h)

print(H)
f = lagr(H)
np.polyval(f, 1)
>>> -20992.0

调用np.polyval(f,1)得到-20992.0，预期结果应为27，需排查问题原因并解决。

问题原因分析 #

• 高次多项式的数值不稳定性（龙格现象）：拉格朗日插值在处理高次多项式（如59阶）时，会出现龙格现象——多项式在插值区间边缘剧烈震荡，导致求值结果严重偏离真实值。尤其当插值点为等间距分布时，这种震荡会被放大。
• 浮点数运算累积误差：传统拉格朗日插值计算基函数时，分母是多个差值的乘积，高次场景下这个乘积的数值会变得极大或极小，浮点数的有限精度无法承载，累积误差迅速放大，最终导致多项式系数完全失真。
• 插值点x值的等间距分布：测试数据中H的x值是1到60的连续整数，属于等间距分布，这种分布在高次拉格朗日插值中是最容易触发龙格现象的情况之一。

解决方案 #

1. 改用数值稳定性更好的重心形式拉格朗日插值 #

重心形式的拉格朗日插值避免了直接构建高次多项式，而是通过计算重心权重来直接求值，数值稳定性远优于传统实现。代码示例：

import numpy as np

def lagrange_barycentric(P):
    x = np.array([p[0] for p in P])
    y = np.array([p[1] for p in P])
    n = len(x)
    # 计算重心权重
    w = np.ones(n)
    for i in range(n):
        # 排除当前点的所有x值
        temp = x[i] - x[np.arange(n) != i]
        w[i] = 1.0 / np.prod(temp)
    
    # 定义求值函数
    def interpolate(x_eval):
        # 避免除以0（当x_eval等于某个插值点x时）
        mask = x_eval != x
        if not np.all(mask):
            return y[np.argwhere(x == x_eval)[0][0]]
        numerator = np.sum(y * w / (x_eval - x))
        denominator = np.sum(w / (x_eval - x))
        return numerator / denominator
    
    return interpolate

使用该函数测试你的数据：

A = [27, 35.1, 26.3, 32.5, 29.1, 32.9, 28.4, 33.5, 26.4, 32.1, 26, 31.4, 25.7, 30.9, 24, 27.8, 25, 27.6, 24.8, 30.5, 23.9, 29.6,
     27.3, 31.8, 26.1, 30.4, 25.7, 29.5, 26.4, 33, 26.4, 34, 26.9, 32.1, 28, 34.9, 26.5, 34.1, 27, 32.7, 27.2, 31.1, 27.6, 33.8,
     27.9, 34.8, 28.4, 34.2, 27.6, 33.3, 28, 32.8, 28.6, 33.1, 27.9, 32.5, 26.3, 29, 25.9, 31]
B = list(range(1, 61))
H = [[B[i], A[i]] for i in range(len(A))]

f = lagrange_barycentric(H)
print(f(1))  # 输出接近27，符合预期

2. 替换为低次分段插值方法 #

如果不需要严格的高次多项式拟合，推荐使用三次样条插值或分段线性插值，这类方法通过将区间拆分为多个小段，每段用低次多项式拟合，完全避免高次插值的数值不稳定问题。示例（使用scipy的三次样条）：

from scipy.interpolate import CubicSpline

x = np.array(B)
y = np.array(A)
cs = CubicSpline(x, y)
print(cs(1))  # 输出27，完全符合预期

3. 调整插值点分布（可选） #

如果必须使用高次插值，将等间距的插值点替换为切比雪夫节点分布，可以显著缓解龙格现象。切比雪夫节点在区间内的分布是两端密集、中间稀疏，能有效降低多项式震荡。

内容的提问来源于stack exchange，提问作者korean-05