3

IDLからPythonにコードを移動する作業をしています。IDL呼び出しの1つは、固定範囲で統合を実行するINT_TABULATEへの呼び出しです。

INT_TABULATED関数は、5点のニュートン・コーツ積分式を使用して、閉区間[MIN(x)、MAX(x)]で表形式のデータセット{xi、fi}を積分します。

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

結果は曲線の下の領域です。

IDLドキュメント

私の質問は、Numpy / SciPyは同様の形式の統合を提供しますか?それは存在するように見えますが、 「面積ではなくニュートン・コーツ積分の重みと誤差係数[scipy.integrate.newton_cotes]」を返すように見えます。

4

1 に答える 1

6

Scipy は、既定では、表形式のデータに対してこのような高次のインテグレーターを提供していません。自分でコーディングせずに利用できる最も近いscipy.integrate.simpsものは、3 点ニュートン コーツ法を使用する です。

単純に同等の積分精度を得たい場合は、xf配列を 5 つのポイントのチャンクに分割しscipy.integrate.newton_cotes、次の行に沿って何かを実行して返される重みを使用して、一度に 1 つずつ積分することができます。

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 点ニュートン コートを実行します。データが等間隔でない場合は、3 次スプラインを使用してそれらの点の値を補間し、等間隔グリッドを構築します。コードはこちらで確認できます。

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
于 2013-01-15T21:35:04.223 に答える