IDL 的 INT_TABULATE - SciPy 等价物？

Question

我正在将一些代码从 IDL 移到 python 中。一个 IDL 调用是 INT_TABULATE，它在固定范围内执行积分。

INT_TABULATED 函数使用五点 Newton-Cotes 积分公式对闭区间 [MIN(x) , MAX(x)] 上的一组表格数据 { xi , fi } 进行积分。

Result = INT_TABULATED( X, F [, /DOUBLE] [, /SORT] )

结果是曲线下的面积。

IDL 文档

我的问题是，Numpy/SciPy 是否提供了类似的集成形式？我看到它[scipy.integrate.newton_cotes]存在，但它似乎返回“牛顿-科特斯积分而不是面积的权重和误差系数”。

Answer 1

默认情况下，Scipy 不为列表数据提供如此高阶的积分器。无需自己编码即可获得的最接近的方法是scipy.integrate.simps，它使用 3 点 Newton-Cotes 方法。

如果您只是想获得可比较的积分精度，您可以将您的x和f数组拆分为 5 个点块，并使用通过scipy.integrate.newton_cotes执行以下操作返回的权重一次整合一个：

def idl_tabulate(x, f, p=5) :
    def newton_cotes(x, f) :
        if x.shape[0] < 2 :
            return 0
        rn = (x.shape[0] - 1) * (x - x[0]) / (x[-1] - x[0])
        weights = scipy.integrate.newton_cotes(rn)[0]
        return (x[-1] - x[0]) / (x.shape[0] - 1) * np.dot(weights, f)
    ret = 0
    for idx in xrange(0, x.shape[0], p - 1) :
        ret += newton_cotes(x[idx:idx + p], f[idx:idx + p])
    return ret

这会在所有时间间隔上执行 5 点牛顿-柯特，除了最后一个，它将执行剩余点数的牛顿-柯特。不幸的是，这不会给你相同的结果，IDL_TABULATE因为内部方法不同：

• Scipy 使用看起来像最小二乘拟合的方法计算不等间距的点的权重，我不完全明白发生了什么，但代码是纯 python，你可以在你的 Scipy 安装文件中找到它scipy\integrate\quadrature.py。

• INT_TABULATED始终对等距数据执行 5 点 Newton-Cotes。如果数据不是等距的，它会构建一个等距的网格，使用三次样条对这些点处的值进行插值。您可以在此处查看代码。

对于INT_TABULATEDdocstring 中的示例，它应该1.6271使用原始代码返回，并且具有 的精确解1.6405，上述函数返回：

>>> x = np.array([0.0, 0.12, 0.22, 0.32, 0.36, 0.40, 0.44, 0.54, 0.64,
...               0.70, 0.80])
>>> f = np.array([0.200000, 1.30973, 1.30524, 1.74339, 2.07490, 2.45600,
...               2.84299, 3.50730, 3.18194, 2.36302, 0.231964])
>>> idl_tabulate(x, f)
1.641998154242472