数値積分

この記事では、SciPyを活用した数値積分の方法について解説します。基本的な一変数積分から始めて、広義積分、重積分、フレネル積分といった応用的な積分についても豊富なサンプルコードを掲載しています。記事後半ではSymPyによるシンボリック積分についても解説します。

【SciPy】数値積分ライブラリ

SciPyの積分パッケージ scipy.integrate.quad() は Fortran のライブラリ QUADPACK を使って、数値積分を実行します。必須引数は被積分関数 func、積分下限値 a、積分上限値 b です。戻り値は積分近似値と推定誤差のタプルです。
 
例として、$y=\sin^2 x$ を $0$ から $\pi$ まで積分してみます (積分の真値は $\pi/2$)。

# PYTHON_INTEGRATE_01

# In[1]

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

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

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

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

print(integ)
# (1.5707963267948966, 1.743934249004316e-14)

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

# In[2]

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

# yを0から∞まで積分して積分値と推定誤差を求める
val, err = integrate.quad(y, 0, np.inf)

print(val)
print(err)

# 1.0000000000000002
# 5.842607038578007e-11

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

# In[3]

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

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

print(val)
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$ です)。
 
方法の１つは integrate.quad() の オプション引数 weight に “cauchy” を渡すことです。weight は被積分関数 func に乗じる重みを表す引数です。weight=”cauchy” は 1/(x-c) を掛けることを意味し、分母の定数 c は wvar で指定します。

# PYTHON_INTEGRATE_02

# In[1]

from scipy import integrate

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

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

print(val)
# 2.220446049250313e-16

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

# In[2]

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

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

print(val)
# 0.0

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

重積分

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$)。

# DOUBLE_INTEGRAL

# In[1]

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

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

# fを0からpiまで数値積分
val, err = integrate.dblquad(f, 0, 1, 0, 1)

print(val)
print(err)

# 0.16666666666666669
# 5.527033708952211e-15

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\]
を計算する場合は次のようなコードを書きます。

# In[2]

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

# fを三重積分
val, err = integrate.tplquad(f, 0, 1, 0, 1, 0, 1)

print(val)
print(err)

# 1.5
# 2.7707360439619496e-14

台形公式

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

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

台形公式は、下図のように積分区間を等分割し、複数の台形を寄せ集めて積分を近似します。
 
 
簡単な連続関数であれば、台形公式で十分な精度の値が得られることが知られています。また、単純なアルゴリズムで処理するので、scipy.integrate.quad() に比べて計算時間が短くて済むという利点もあります。
 
一例として、二次関数 $y=x^2$ を $0$ から $1$ まで積分してみます。

# INTEGRATE_TRAPEZOIDAL

# In[1]

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() は シンプソンの公式 (Simpson’s rule) にしたがって数値積分を実行します。

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

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

# INTEGRATE_SIMPSON

# In[1]

# 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

フレネル積分

フレネル正弦積分とフレネル余弦積分は
 \[\mathrm{S}\,(x)=\int_{0}^{x}\sin\left(\frac{\pi}{2}t^2\right)dt,\quad\mathrm{C}\,(x)=\int_{0}^{x}\cos\left(\frac{\pi}{2}t^2\right)dt\]
によって定義され、まとめてフレネル積分とよばれます。フレネル積分は初等関数で表すことができません（級数展開で表記することは可能です）。scipy.special.fresnel(z) は z におけるフレネル正弦積分とフレネル余弦積分の値を返します。以下のコードを実行するとフレネル積分の概形を描きます。

# FRESNEL_INTEGRAL

# In[1]

import numpy as np
from scipy.special import fresnel
import matplotlib.pyplot as plt

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

# フレネル正弦積分
fs = fresnel(x)[0]

# フレネル余弦積分
fc = fresnel(x)[1]

# FigureとAxes
fig = plt.figure(figsize=(6.5, 5))
ax = fig.add_subplot(111)
ax.grid()
ax.set_xlabel("x", fontsize=15)
ax.set_ylabel("S(x), C(x)", fontsize=15)
ax.set_xlim(-4, 4)
ax.set_ylim(-1, 1)

# Axesにグラフをプロット
ax.plot(x, fs, color = "blue", label="S (x)")
ax.plot(x, fc, color = "red", label="C (x)")

# 凡例を表示
ax.legend(loc="upper left")

plt.show()

指数積分

scipy.special.expi(x) は指数積分
 \[\mathrm{Ei}(x)=\int_{-\infty}^{x}\frac{e^{-t}}{t}dt\]
の値を返します。以下のコードを実行すると、実変数の $\mathrm{Ei}(x)$ のグラフの概形を描きます。

# EI_INTEGRAL

# In[1]

import numpy as np
from scipy.special import expi
import matplotlib.pyplot as plt

# 区間[-4,4]を1024分割
x = np.linspace(-4, 4, 1025)

# 指数積分Ei(x)
ei = expi(x)

# FigureとAxes
fig = plt.figure(figsize=(6.5, 5))
ax = fig.add_subplot(111)
ax.grid()
ax.set_xlabel("x", fontsize=15)
ax.set_ylabel("Ei(x)", fontsize=15)
ax.set_xlim(-4, 4)
ax.set_ylim(-4, 8)

# 指数積分をプロット
ax.plot(x, ei, color = "blue", label="Ei(x)")

# 凡例を表示
ax.legend(loc="upper left")

plt.show()

 
scipy.special.expn(x) は $n$ 次の指数積分
 
\[\int_{1}^{\infty}\frac{e^{-xt}}{t^n}dt\]
の値を返します。

# In[2]

from scipy.special import expn

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

# FigureとAxes
fig = plt.figure(figsize=(6.5, 5))
ax = fig.add_subplot(111)
ax.grid()
ax.set_xlabel("x", fontsize=15)
ax.set_ylabel("En(x)", fontsize=15)
ax.set_xlim(0, 2)
ax.set_ylim(0, 10)

# 指数積分をプロット
for n in range(3):
    ax.plot(x, expn(n, x), label = "n={}".format(n))

# 凡例を表示
ax.legend(loc="upper right")

plt.show()

 
scipy.special.exp1(x,n) は指数積分
 \[\int_{1}^{\infty}\frac{e^{-xt}}{t}dt\]
の値を返します。scipy.special.expn(x,n) で n = 1 を指定した場合の戻り値と同じです。

ドーソン積分

scipy.special.dawsn(x) はドーソン積分 (Dawson integral)
 \[D_+(x)=\exp(-x^2)\int_0^x\exp(t^2)dt\]
を計算します。

# DAWSON_INTEGRAL

# In[1]

import numpy as np
from scipy.special import dawsn
import matplotlib.pyplot as plt

# 区間[-10,10]を512分割
x = np.linspace(-10, 10, 513)

# 正弦積分
ds = dawsn(x)

# FigureとAxes
fig = plt.figure(figsize=(6.5, 5))
ax = fig.add_subplot(111)
ax.grid()
ax.set_title("Dawson integral", fontsize=16)
ax.set_xlabel("x", fontsize=15)
ax.set_ylabel("D(x)", fontsize=15)
ax.set_xlim(-10, 10)
ax.set_ylim(-0.6, 0.6)

# グラフをプロット
ax.plot(x, ds, color = "blue")

plt.show()

【SymPy】不定積分と定積分

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

integrate(func, var)

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

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

例として、対数関数 $\log x$ の不定積分公式
 \[\int \log{x}dx = x\log x-x+C\]
を確認してみましょう。

# SYMPY_INTEGRATE

# In[1]

import sympy

# 数式をLaTeXで表示
sympy.init_printing()

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

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

# 対数関数の不定積分
integ = sympy.integrate(fx, x)

display(integ)

sympy.integrate() の上限と下限を指定すると定積分を計算します。$\log x$ の $x=1$ から $x=2$ までの定積分を求めてみましょう。

# In[2]

# logxを1から2まで定積分
integ = sympy.integrate(fx, (x, 1, 2))

display(integ)

そのままでは、$\log(2)$ などは計算しないので、値がほしいときには evalf() メソッドで評価値を得てください。

# In[3]

# 定積分の評価値を取得
val = integ.evalf(5)

print(val)
# 0.38629

上限もしくは下限に無限量を指定したいときには、sympy.oo を使います。たとえば、
 \[\int_{-\infty}^{\infty}\frac{1}{1+x^2}dx\]
を計算する場合は以下のようなコードを書きます。

# In[4]

# 1/(1+x**2)を-∞から+∞まで定積分
integ = sympy.integrate(1/(1+x**2), (x, -sympy.oo, sympy.oo))

display(integ)

sympy.integrate() を使うと、手計算では扱いにくいような複雑な積分を簡単に実行できます。一例として、$\tan^3 x$ を積分してみましょう。

# In[5]

# tanxの3乗の不定積分
integ = sympy.integrate(sympy.tan(x)**3)
integ = sympy.simplify(integ)

display(integ)