INT_TABULATE - SciPy的等价物？

Question

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

INT_TABULATED函数在闭区间MIN(x)，MAX(x)上集成了一组表列数据{ xi，fi }，使用五点Newton-Cotes积分公式。

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

结果是曲线下的面积。

IDL DOCS

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

Answer 1

在默认情况下，Scipy不为列表数据提供这样一个高阶积分器。最近的你没有编码它自己可用的是scipy.integrate.simps，它使用3点牛顿-科特斯方法。

如果您只是想获得类似的集成精度，可以将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相同的结果，因为内部方法是不同的：

• 使用最小面积的合适值来计算不相等间距的点的权重，我不完全理解发生了什么，但是代码是纯python，您可以在文件scipy\integrate\quadrature.py中找到它。
• INT_TABULATED总是对等速数据执行5点牛顿-科特函数.如果数据不是等距的，则使用三次样条插值这些点的值，建立等距网格。您可以检查代码这里。

对于INT_TABULATED docstring中的示例，它可以使用原始代码返回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