SciPyによる積分(scipy.integrate)①｜SciPy入門 #8
SciPyについて色々と話題になり面白そうだったので公式チュートリアルを元にまとめています。
SciPy Tutorial — SciPy v1.2.1 Reference Guide
#5、#6、#7ではscipy.optimizeを元に最適化に関して取り扱いました。
#7でscipy.statsのチュートリアルは一通り終了したので、#8からは積分に関する機能であるscipy.integrateについて取り扱っていきます。
Integration (scipy.integrate) — SciPy v1.2.1 Reference Guide
#8ではチュートリアルからGeneral integrationとintegrationの二つについて取り扱います。以下目次になります。
0. モジュール概要
1. General integration (quad)
2. integration (dblquad, tplquad, nquad)
3. まとめ
0. モジュール概要
モジュール概要に関してはhelp(integrate)を実行することで確認することができます。
上記のように、quadに始まり様々な機能が実装されています。
1. General integration (quad)
まずscipy.integrate.quadについて取り扱っていきます。こちらのメソッドを使用することで、1変数関数の2点間での積分を行うことができます。
import scipy.integrate as integrate
result = integrate.quad(lambda x: x, 0, 4)
print(result[0])
上記の実行例ではを区間[0, 4]で積分したものが結果となっています。短い辺の長さが4の直角二等辺三角形の面積に一致するので、この値は
となります。上記の実行例では積分結果のresult[0]として8が出力されるようになっています。
次に積分実行時に関数にパラメータを与えるパターンの解説がされています。
from scipy.integrate import quad
def integrand(x, a, b):
return a*x + ba = 2
for b in range(0,5):
I = quad(integrand, 0, 1, args=(a,b))
print(str(b)+": "+str(I[0]))
上記はチュートリアルのコードを元にわかりやすくするために少々変更を行いました。傾きが2の一次関数において、切片を0~4に変更して、積分値を計算しています。
実行結果は上記のように1,2,3,4,5と切片が1増加するにつれて積分値も1ずつ増加しており、正しい結果であるということがわかります。
チュートリアルでは上記の数式の積分を計算しています。実行コードは以下になります。
from scipy.integrate import quad
def integrand(t, n, x):
return np.exp(-x*t) / t**n
def expint(n, x):
return quad(integrand, 1, np.inf, args=(n, x))[0]
result = quad(lambda x: expint(3, x), 0, np.inf)
print(result[0])
scipy.special.expnでdefで定義したexpintと同様の関数は実装されているのですが、quadの引数として与えるにあたって別途定義をしている形になります。積分区間は[0,∞)を引数として与えており、無限大はnp.infで与えていることに注目です。実行結果は下記のようになります。
途中関数expintの値の検証を行うにあたって、np.vectorizeを用いてscipy.special.expnの値と比較を行っていますが結果には影響ありません。結果は0.3333..となり、これが積分を行った結果となります。
2. integration (dblquad, tplquad, nquad)
重積分について諸々まとまっています。
使用例として上記の式の積分についてがまとまっています。
from scipy.integrate import quad, dblquad
def I(n):
return dblquad(lambda t, x: np.exp(-x*t)/t**n, 0, np.inf, lambda x: 1, lambda x: np.inf)
print(I(4)[0])
print(I(3)[0])
print(I(2)[0])
上記を実行すると0.2500.., 0.333.., 0.4999..という結果が得られ、数式と同様の結果が得られていることがわかります。また、ここで用いたdblquadは二重積分を行うための関数です。この辺は少し難易度が上がるので今回は一例に止めたいと思います。
3. まとめ
#8ではscipy.integrate.quadやscipy.integrate.dblquadなどについて取り扱いました。
#9ではscipy.integrateのチュートリアルの残りについてまとめたいと思います。
