python - IDL 的 INT_TABULATE - SciPy 等价物？

Question

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

The INT_TABULATED function integrates a tabulated set of data { xi , fi } on the closed interval [MIN(x) , MAX(x)], using a five-point Newton-Cotes integration formula.

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

其中结果是曲线下的面积。

IDL DOCS

我的问题是，Numpy/SciPy 是否提供类似形式的集成？我看到 [scipy.integrate.newton_cotes] 存在，但它似乎返回“weights and error coefficient for Newton-Cotes integration 而不是区域”。

Answer 1

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

如果您只是想获得可比较的积分精度，您可以将 x 和 f 数组拆分为 5 个点 block ，并使用权重一次对它们进行积分由 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 点 Newton-Cotes，也许除了最后一个，它将对剩余点数执行 Newton-Cotes。不幸的是，这不会给您与 IDL_TABULATE 相同的结果，因为内部方法不同:

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

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

对于 INT_TABULATED 文档字符串中的示例，它应该使用原始代码返回 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