3

我正在将一些代码从 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]存在，但它似乎返回“ Newton-Cotes 积分而不是面积的权重和误差系数”。

4

1 回答 1

6

默认情况下，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 点 Newton-Cotes，除了最后一个，它将对剩余的点数进行 Newton-Cotes。不幸的是，这不会给您相同的结果，IDL_TABULATE因为内部方法不同：

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

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

对于文档字符串中的示例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
于 2013-01-15T21:35:04.223 回答