SciPy 数値積分

SciPy 数値積分

SciPy 数値積分

scipy.integrate.quad()

 scipy.integrate.quad() は Fortran のライブラリ QUADPACK を使って、ガウスの数値積分 (Gaussian quadrature) を実行します。必須引数は被積分関数 func、積分下限値 a、積分上限値 b です。戻り値は積分近似値と推定誤差のタプルです。

 例として、$y=\sin^2 x$ を $0$ から $\pi$ まで積分してみます (積分の真値は $\pi/2$)。

# NI_01-1

# NumPyをインポート
import numpy as np

# SciPy積分パッケージをインポート
from scipy import integrate

# 関数を定義
y = lambda x: np.sin(x)**2

# yを0からpiまで数値積分
ie = integrate.quad(y, 0, np.pi)

print(ie)
(1.5707963267948966, 1.743934249004316e-14)

 

無限区間における積分

 下限 a または上限 b に $\infty$ を指定することもできます。
 たとえば、$\displaystyle\int_0^{\infty}e^{-x}dx$ を実行する場合は以下のように記述します (真値は $1$)。

# NI_01-2

# 関数を定義
y = lambda x: np.exp(-x)

# yを0から∞まで積分して積分値と推定誤差をiyとerrに格納
iy, err = integrate.quad(y, 0, np.inf)

print(iy)
print(err)
1.0000000000000002
5.842607038578007e-11

 

パラメータの設定

 関数が積分変数以外のパラメータを含む場合は、args でパラメータの値を設定できます。たとえば、
 
\[\int_{0}^{1}\frac{dx}{x^2+a^2}\quad (a=2)\]
を計算する場合は以下のように記述します。

# NI_01-3

# 関数を定義
y = lambda x, a: 1/(x**2 + a**2)

# a=2として、yを0から1まで積分
iy, err = integrate.quad(y, 0, 1, args=2)

print(iy)
print(err)
0.23182380450040307
2.5737612541254846e-15

 

特異点を含む関数の積分

 分数関数 $f(x)=1/x$ の区間 $[-1,1]$ における積分
 
\[\int_{-1}^{1}\frac{1}{x}dx\]
は $x=0$ で $f(x)$ が特異点となるので、通常の積分は値をもちませんが、特異点を除いた積分 (コーシーの主値積分)
 
\[\lim_{\varepsilon\rightarrow +0}\left\{\int_{a}^{c-\varepsilon}f(x)dx+\int_{c+\varepsilon}^{b}f(x)dx\right\}\]
を計算することは可能です ($1/x$ は原点に関して対称なので積分の真値は $0$ です)。

 方法の1つは integrate.quad() の オプション引数 weight に "cauchy" を渡すことです。weight は被積分関数 func に乗じる重みを表す引数です。weight="cauchy" は 1/(x-c) を掛けることを意味し、分母の定数 c は wvar で指定します。

# 関数を定義
y = lambda x : 1

# 1/xのコーシー主値積分
iy, err = integrate.quad(y, -1, 1,
                         weight="cauchy", wvar=0)

print(iy)
2.220446049250313e-16

 もう1つの方法は points に特異点を渡すことです。

# 関数を定義
y = lambda x : 1/x

# 1/xのコーシー主値積分
iy, err = integrate.quad(y, -1, 1, points=0)

print(iy)
0.0

 points には複数の値を渡せるので、積分範囲が複数の特異点を含む場合は、特異点のリストを渡してください。
 

二重積分:scipy.integrate.dblquad()

 scipy.integrate.dblquad() は、二変数関数を受け取って二重積分 (double integral) を実行します。たとえば、領域 $0\leq x,\ y\leq 1$ における $f(x,y)=x^2y$ の積分
 
\[\int_{0}^{1}\int_{0}^{1}x^2ydxdy\]
を計算する場合は以下のコードを実行します (真値は $1/6$)。

# NI_02

# SciPy積分パッケージをインポート
from scipy import integrate

# 関数を定義
f = lambda y, x: (x**2)*y

# yを0からpiまで数値積分
i_f, err = integrate.dblquad(f, 0, 1, lambda x: 0, lambda x: 1)

print(i_f)
print(err)
0.16666666666666669
5.527033708952211e-15

 

三重積分:scipy.integrate.tplquad()

 scipy.integrate.tplquad() に三変数関数を与えると、三重積分 (triple integral) を計算できます。たとえば、$f(x,y,z)=x+y+z$ の空間領域 $0\leq x,\ y,\ z \leq 1$ における積分
 
\[\int_{0}^{1}\int_{0}^{1}\int_{0}^{1}(x+y+z)dxdydz\]
を計算する場合は次のようなコードを書きます。

# NI_03

# SciPy積分パッケージをインポート
from scipy import integrate

# 関数を定義
f = lambda z, y, x: x + y + z

# fを三重積分
i_f, err = integrate.tplquad(f, 0, 1,
                             lambda x: 0, lambda x: 1,
                             lambda x, y: 0, lambda x, y: 1)

print(i_f)
print(err)
1.5
2.7707360439619496e-14

 

入門 Python 3

新品価格
¥3,996から
(2019/8/6 11:45時点)

台形公式:scipy.integrate.trapz()

 scipy.integrate.trapz() は 台形公式 (trapezoidal rule) を使って積分の近似値を求めます。

 scipy.integrate.trapz(y, x=None, dx=1.0, axis=-1)

 台形公式は、下図のように積分区間を等分割し、複数の台形を寄せ集めて積分を近似します。

 Scipy 台形公式 (trapezoidal rule) による数値積分

 簡単な連続関数であれば、台形公式で十分な精度の値が得られることが知られています。
 また、単純なアルゴリズムで処理するので、scipy.integrate.quad() に比べて計算時間が短くて済むという利点もあります。

 一例として、二次関数 $y=x^2$ を $0$ から $1$ まで積分してみます。

# NI_04

import numpy as np

# 積分ライブラリをインポート
from scipy import integrate

# 区間[0,1]を64分割
x = np.linspace(0, 1, 65)

# xの各要素に対するyを計算
y = x**2

# 台形公式を使ってyを0から1まで積分
s = integrate.trapz(y, x)

print(s)
0.3333740234375

 区間を $64$ 分割して計算させています。真値は $1/3=0.33333...$ なので若干の誤差がありますが、x の刻み幅を小さくすることで精度を向上させることができます。
 

シンプソンの公式:scipy.integrate.simps()

 scipy.integrate.simps() は シンプソンの公式 (Simpson's rule) にしたがって数値積分を実行します。

 scipy.integrate.simps(y, x=None, dx=1, axis=-1, even='avg')

 シンプソンの公式は曲線を 3 点を通る放物線の集合によって近似します。
 以下のコードでは $y=\sqrt{x}$ を $0$ から $1$ まで積分します。

# NI_05

# NumPyをインポート
import numpy as np

# 積分パッケージをインポート
from scipy import integrate

# 区間[0,1]を128分割
x = np.linspace(0, 1, 129)

# xの各要素に対するyを計算
y = np.sqrt(x)

# シンプソンの公式を使ってyを0から1まで積分
s = integrate.simps(y, x)

print(s)
0.6666106059362655

 

SymPy 積分

 SymPy では次の構文で関数 func の変数 var による不定積分を求めることができます。

 integrate(func, var)

 定積分を計算する場合は次のように記述します。

 integrate(func, (var, 下限, 上限))

 次のサンプルコードでは、$\log x$ の不定積分と定積分を計算しています。

# SI_01

import sympy

# 記号a, b, xを定義
sympy.var('a b x')

# f(x) = log(x)
f = sympy.log(x)

# logxの不定積分
i = sympy.integrate(f, x)

# logxをaからbまで積分
i_a_b = sympy.integrate(f, (x, a, b))

# logxを1から2まで積分
i_1_2 = sympy.integrate(f, (x, 1, 2))

print("logxの不定積分 {}".format(i))
print("logxをaからbまで積分 {}".format(i_a_b))
print("logxを1から2まで積分 {}".format(i_1_2))
logxの不定積分 x*log(x) - x
logxをaからbまで積分 -a*log(a) + a + b*log(b) - b
logxを1から2まで積分 -1 + 2*log(2)

 変数 var の下限と上限に oo のような特殊な記号を指定して広義積分を求めることもできます。