Python Numpy与Excel六阶多项式计算结果不一致的原因排查

运行环境：Windows 11 x64系统

背景 #

在Excel中对数据做六阶多项式线性回归，得到如下系数：

a1 =  -0.000000000000000000051066485848517600000
a2 = 0.0000000000000012806568632620900000
a3 = -0.00000000000531035856252443000
a4 = 0.0000000086330792556790300
a5 = -0.0000063412693200472300000
a6 = 0.0020307237682114700000000
b1 = 19.2031257127800000000000000

对应的多项式公式为：
a1*f**6 + a2*f**5 + a3*f**4 + a4*f**3 + a5*f**2 + a6*f + b1

在50至1800的f取值范围内，Excel计算出的曲线符合预期（趋势平缓，随f变化呈现合理渐变），但Python运行下方代码得到的曲线结果与Excel完全不符（出现异常波动或趋势偏离）。

Python代码 #

import numpy as np
import matplotlib.pylab as plt

f = np.arange(50,1801)

def gain_response (f):
    a1 =  np.longdouble(-0.000000000000000000051066485848517600000)
    a2 = np.longdouble(0.0000000000000012806568632620900000)
    a3 = np.longdouble(-0.00000000000531035856252443000)
    a4 =  np.longdouble(0.0000000086330792556790300)
    a5 =  np.longdouble(-0.0000063412693200472300000)
    a6 = np.longdouble(0.0020307237682114700000000)
    b1 = np.longdouble(19.2031257127800000000000000)
    
    return a1*f**6 + a2*f**5 + a3*f**4 + a4*f**3 + a5*f**2 + a6*f + b1
    
plt.plot(f, gain_response(f), label="-40C")
plt.legend()
plt.grid()
plt.tight_layout()
plt.show()

问题原因 #

核心问题不是np.longdouble的精度不够，而是高次多项式直接展开求值的数值稳定性极差：

• 当f取到1800时，f**6是约3.4e22的超大数，乘以极小的a1（约-5e-20）后，计算过程会损失大量有效数字；
• 各项量级差异极大（比如b1是19，a6*f是3600左右），直接相加时，小量级的项会被大量级的项“吞噬”，浮点数舍入误差被放大，导致结果完全偏离预期。

而Excel默认使用**霍纳法则（Horner's Method）**优化多项式求值，通过嵌套乘法加法的形式，极大提升了数值稳定性，减少了精度损失。

解决方案 #

在Python中改用霍纳法则改写计算逻辑，将原式转换为嵌套形式：
((((a1*f + a2)*f + a3)*f + a4)*f + a5)*f + a6)*f + b1

修改后的代码如下：

import numpy as np
import matplotlib.pylab as plt

f = np.arange(50,1801)

def gain_response (f):
    a1 =  np.longdouble(-0.000000000000000000051066485848517600000)
    a2 = np.longdouble(0.0000000000000012806568632620900000)
    a3 = np.longdouble(-0.00000000000531035856252443000)
    a4 =  np.longdouble(0.0000000086330792556790300)
    a5 =  np.longdouble(-0.0000063412693200472300000)
    a6 = np.longdouble(0.0020307237682114700000000)
    b1 = np.longdouble(19.2031257127800000000000000)
    
    # 霍纳法则求值
    return (((((a1 * f + a2) * f + a3) * f + a4) * f + a5) * f + a6) * f + b1
    
plt.plot(f, gain_response(f), label="-40C")
plt.legend()
plt.grid()
plt.tight_layout()
plt.show()

修改后重新运行，计算结果会和Excel一致，曲线趋势恢复正常。

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