Scipy は、既定では、表形式のデータに対してこのような高次のインテグレーターを提供していません。自分でコーディングせずに利用できる最も近いscipy.integrate.simps
ものは、3 点ニュートン コーツ法を使用する です。
単純に同等の積分精度を得たい場合は、x
とf
配列を 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_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