显式传入采样点x时scipy.integrate.simpson结果异常原因问询

问题背景 #

我基于辛普森法则实现了自定义积分函数，将结果与解析解及scipy.integrate.simpson的计算结果对比：

• 在步长为1、区间[0,100]（101个采样点）积分x²时，显式传入x数组的scipy结果为333333.3333333334，与解析解、仅传y的scipy结果、自定义实现结果（均为333333.3333333333）不符；
• 积分3*x²时所有结果一致。

已知该函数支持非偶数区间与非均匀采样点，推测差异可能与此相关，求显式传x导致差异的原因。

测试代码 #

from scipy.integrate import simpson
import numpy as np
x = np.linspace(0,100,101)

y = x**2
y_int = (x**3)/3

y_analytic = np.max(y_int) - np.min(y_int)

sci_int = simpson(y,x) #sampled points explicitly given
sci_int_2 = simpson(y) #assumes dx=1

def custom_simpson(x,y):
    a = x[0]
    b = x[-1]

    n = len(x) - 1 #number of intervals

    h = (b-a) / n

    even_y = y[::2]
    odd_y = y[1::2]

    even_y = even_y[1:-1] 
    
    odd_y = np.sum(odd_y)
    even_y = np.sum(even_y)
    
    result = (h/3) * (y[0] + y[-1] + (4*odd_y) + (2*even_y))

    return result


print('hand_made \t', custom_simpson(x,y))
print('analytic \t', y_analytic)
print('scipy_dx=1 \t', sci_int_2)
print('scipy_x_given \t', sci_int)

测试输出 #

当y = x**2时 #

hand_made        333333.3333333333
analytic         333333.3333333333
scipy_dx=1       333333.3333333333
scipy_x_given    333333.3333333334

当y = 3*(x**2)时 #

hand_made        1000000.0
analytic         1000000.0
scipy_dx=1       1000000.0
scipy_x_given    1000000.0

原因分析 #

这个差异源于浮点数计算的精度误差，以及scipy.integrate.simpson在不同参数传入方式下的计算路径差异：

• 当显式传入x数组时，scipy会基于每个相邻x点的间隔计算积分（即使间隔均匀）；而默认传入y时，会直接使用固定步长dx=1计算。对于均匀采样的x数组，两种方式数学上等价，但浮点数运算的累积误差会导致细微差异。
• 积分x²的解析解为100³/3 ≈ 333333.3333333333，这个值无法用二进制浮点数精确表示，会存在舍入误差。scipy在显式处理x数组时，内部数值计算步骤（如逐个间隔处理）导致舍入后的结果与其他实现略有不同，表现为末尾的3334。
• 积分3*x²的解析解为1000000.0，这个值是2的整数次幂的倍数，可以被二进制浮点数精确表示，因此所有计算路径的结果完全一致，不会出现精度差异。

补充：你的自定义实现和scipy默认dx=1的计算逻辑更接近，都是直接使用整体步长h=(b-a)/n；而显式传x时，scipy的内部实现会遍历每个相邻间隔并应用辛普森公式，虽然数学等价，但浮点数运算的步骤差异导致了最终结果的微小误差。

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