scipy.integrate との連携｜Python
SciPy の integrate モジュールは、数値積分のための強力な関数群を提供する。自前で実装するより高精度で高速であり、実務では SciPy を使うのが一般的だ。主要な関数の使い方を押さえておこう。
quad：汎用的な定積分
quad は最もよく使う関数で、適応型のガウス求積法を内部で使用する。戻り値は積分値と誤差の推定値のタプルだ。
from scipy.integrate import quad
import math
# ∫₀^π sin(x) dx = 2
result, error = quad(lambda x: math.sin(x), 0, math.pi)
print(f'積分値: {result}') # 2.0
print(f'誤差推定: {error}') # 2.2e-14無限区間の積分
quad は無限区間の積分もサポートする。np.inf または math.inf を境界に指定する。
from scipy.integrate import quad
import numpy as np
# ∫₀^∞ e^(-x²) dx = √π / 2
result, error = quad(lambda x: np.exp(-x**2), 0, np.inf)
print(f'積分値: {result}') # 0.8862269254527579
print(f'√π / 2: {np.sqrt(np.pi)/2}') # 0.8862269254527579負の無限大も同様に扱える。
from scipy.integrate import quad
import numpy as np
# ∫_{-∞}^{∞} e^(-x²) dx = √π
result, error = quad(lambda x: np.exp(-x**2), -np.inf, np.inf)
print(f'積分値: {result}') # 1.7724538509055159特異点を含む積分
被積分関数が境界で特異点を持つ場合、quad は自動的に対処しようとする。ただし、特異点の位置を points 引数で明示すると精度が向上することがある。
from scipy.integrate import quad
import numpy as np
# ∫₀^1 1/√x dx = 2（x=0 で発散）
result, error = quad(lambda x: 1/np.sqrt(x), 0, 1)
print(f'積分値: {result}') # 1.9999999999999998dblquad：二重積分
二次元領域での積分には dblquad を使う。内側の積分範囲は外側の変数の関数として指定できる。
from scipy.integrate import dblquad
# ∫₀^1 ∫₀^1 x*y dy dx = 1/4
result, error = dblquad(
lambda y, x: x * y, # 注意: 引数の順序は (y, x)
0, 1, # x の範囲
0, 1 # y の範囲（定数の場合はそのまま数値を渡す）
)
print(f'積分値: {result}') # 0.25引数の順序が (y, x) である点に注意。これは内側の積分（）を先に行うためだ。
tplquad：三重積分
三次元領域での積分には tplquad を使う。
from scipy.integrate import tplquad
# 単位立方体上の積分 ∫∫∫ x*y*z dz dy dx = 1/8
result, error = tplquad(
lambda z, y, x: x * y * z, # 引数順序: (z, y, x)
0, 1, # x の範囲
0, 1, # y の範囲
0, 1 # z の範囲
)
print(f'積分値: {result}') # 0.125nquad：任意次元の積分
4 次元以上の積分や、より柔軟な指定が必要な場合は nquad を使う。
from scipy.integrate import nquad
# 4 次元積分の例
def f(x1, x2, x3, x4):
return x1 * x2 * x3 * x4
ranges = [(0, 1), (0, 1), (0, 1), (0, 1)]
result, error = nquad(f, ranges)
print(f'積分値: {result}') # 0.0625（= 1/16）fixed_quad：固定点ガウス求積法
関数が滑らかで、評価回数を固定したい場合は fixed_quad が効率的だ。
from scipy.integrate import fixed_quad
import numpy as np
# 5 点ガウス・ルジャンドル求積法
result, _ = fixed_quad(np.sin, 0, np.pi, n=5)
print(f'積分値: {result}') # 2.0000000000n 点のガウス求積法は 次までの多項式を正確に積分する。
romberg：ロンバーグ積分
台形則を繰り返し改良するロンバーグ積分も用意されている。
from scipy.integrate import romberg
import numpy as np
result = romberg(np.sin, 0, np.pi, show=True)
# ステップごとの収束の様子が表示されるshow=True で収束過程を確認できる。教育目的に便利だ。
trapezoid と simpson：データ点からの積分
関数ではなくデータ点が与えられている場合は、trapezoid（旧 trapz）や simpson を使う。
from scipy.integrate import trapezoid, simpson
import numpy as np
# データ点
x = np.linspace(0, np.pi, 11)
y = np.sin(x)
# 台形則
print(trapezoid(y, x)) # 1.9835235375094546
# シンプソン則
print(simpson(y, x=x)) # 2.0000067844418012関数を渡す。適応型で高精度。
データ点を渡す。実験データや離散値の積分向き。
solve_ivp：常微分方程式の初期値問題
積分とは少し異なるが、integrate モジュールには常微分方程式のソルバーもある。
from scipy.integrate import solve_ivp
import numpy as np
# dy/dt = -2y, y(0) = 1 の解は y = e^(-2t)
def dydt(t, y):
return -2 * y
sol = solve_ivp(dydt, [0, 3], [1], t_eval=np.linspace(0, 3, 10))
print(sol.y[0]) # 各時刻での y の値odeint（旧 API）より solve_ivp の方が推奨されている。
関数の選び方
一次元の定積分。最初に試すべき関数。
多次元積分。次元が高いと計算コストが急増する。
離散データの積分。実験値や数表からの積分に。
微分方程式の数値解法。動的システムのシミュレーションに。
多くの場合、quad で事足りる。まずは quad を試し、必要に応じて他の関数を検討するとよい。
