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.9999999999999998

dblquad：二重積分

二次元領域での積分には dblquad を使う。内側の積分範囲は外側の変数の関数として指定できる。

$\int_{a}^{b} \int_{g (x)}^{h (x)} f (x, y) d y d x$

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) である点に注意。これは内側の積分（ $y$ ）を先に行うためだ。

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.125

nquad：任意次元の積分

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.0000000000

n 点のガウス求積法は $2 n - 1$ 次までの多項式を正確に積分する。

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

quad 系

関数を渡す。適応型で高精度。

trapezoid / simpson

データ点を渡す。実験データや離散値の積分向き。

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

一次元の定積分。最初に試すべき関数。

dblquad / tplquad / nquad

多次元積分。次元が高いと計算コストが急増する。

trapezoid / simpson

離散データの積分。実験値や数表からの積分に。

solve_ivp

微分方程式の数値解法。動的システムのシミュレーションに。

多くの場合、quad で事足りる。まずは quad を試し、必要に応じて他の関数を検討するとよい。