Python SciPy : SciPy の積分関数の基本的使い方

備忘録として SciPy.integrate パッケージ について、関数オブジェクトの積分関数の基本的使い方をまとめます。

関数の種類

関数オブジェクトの積分関数には以下のようなものがあります。

表1: 関数オブジェクトの積分関数
関数 説明
quad 1 次元積分(適応積分)
dblquad 2 次元積分
tplquad 3 次元積分
nquad 多次元積分
fixed_quad 次数を指定して Gauss 積分
quadrature 絶対誤差の閾値を指定して Gauss 積分
romberg Romberg 積分

1 次元積分をやるときは、とりあえず quad を試せばいいと思います。 これは FORTRAN の積分計算ライブラリ QUADPACK の QAGS と QAGI ルーチンで積分値を計算します。 適応 Gauss-Kronrod 積分なので積分値を比較的少ない演算で精度良く求めることがきます。

2 次元以上の積分は dblquad, tplquad, nquad で計算できます。 これらは内部で 1 次元積分ルーチン(quad)を繰り返し呼び出して計算しています。 適応 Gauss 積分で次数や、絶対誤差の閾値を設定して積分計算をしたい場合はそれぞれ fixed_quad、quadrature で計算します。

ニュートンコーツ系の公式では適応 Romberg 積分が romberg で計算できます。 適応自動積分に適用する場合によく使われるそうです。

使い方の例

関数の使用例をいくつか載せておきます。 どの関数も基本的に被積分関数を定義し、引数に被積分関数と積分区間を指定するだけです。

積分計算で以下のように円周率を求めてみます。

\begin{equation} \pi=\int_0^1\frac{4}{1+x^2}{\rm d}x \nonumber \end{equation}

quad を使って以下のように積分値を求められます。

import numpy as np
from scipy import integrate


def computePi(x):
    return 4/(1+x**2)

integrate.quad(computePi, 0, 1)
(3.1415926535897936, 3.4878684980086326e-14)

結果にはデフォルトで積分値と推定誤差が出力されます。 被積分関数が簡単なら無名関数を用いるとすっきり書くことができます。

integrate.quad(lambda x: 4/(1+x**2), 0, 1)
(3.1415926535897936, 3.4878684980086326e-14)

romberg メソッドも使い方は同じです。

integrate.romberg(computePi, 0, 1)
3.14159265364

積分区間が \([0,\infty]\) の積分は以下のようにして計算できます。

integrate.quad(lambda x: x**5*np.exp(-x)*np.sin(x), 0, np.inf)
(-14.999999999999996, 1.6360613097566042e-08)

関数にパラメータを持たせたい場合は、引数 args に値を設定します。

integrate.quad(lambda x, c : 1/(x**3-2*x-c), 0, 2, args=5)
(-0.4605015338467388, 1.1968913601400428e-08)

更に、いくつかパラメータの値を変更し積分値をベクトルとして求めたい場合は、vectorize メソッドを用いて計算できます。 この辺りは Matlab の integral 関数はもっと簡単に書けるのがいいですね。

def fun(a):
    return integrate.quad(lambda x, a: np.sin(a*x), 0, 1, args=a)[0]

vec_fun = np.vectorize(fun)

a = [1, 2, 3]
vec_fun(a)
array([ 0.45969769,  0.70807342,  0.66333083])

2 次元積分の dblquad も 1 次元積分と同様に使いますが、2 次元目の積分区間は関数で定義するのが注意する所です。 これは 3 次元以上の場合も同様です。 $0≤ x ≤ 1$、\(0\le y \le 1-x\) の三角領域で定義される関数の積分値は以下のように求まります。

integrate.dblquad(lambda x, y: 1/(np.sqrt(x+y)*(1+x+y)**2), 0, 1, lambda x: 0, lambda x: 1-x)
(0.285398163397451, 1.5065093369366904e-08)

コメント

Comments powered by Disqus
書籍更新情報
2017-02-18
Pythonによる科学技術計算 基礎編
1.4版への更新が可能になりました。
サポートページはこちら
電子書籍
Pythonによる科学技術計算 基礎編
電子書籍
線形代数(1): Pythonによる科学技術計算 実践編