数値積分

数値積分

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

 

アルゴリズム図鑑 絵で見てわかる26のアルゴリズム

中古価格
¥2,011から
(2019/7/30 20:26時点)

台形公式: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

 

フレネル積分

 フレネル正弦積分とフレネル余弦積分は
 
\[\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 におけるフレネル正弦積分とフレネル余弦積分の値を返します。以下のコードを実行するとフレネル積分の概形を描きます。

# 区間[-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.savefig("fresnel.png", bbox_inches = "tight")

 Scipy フレネル積分のグラフ
 

指数積分

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

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.savefig("expi.png", bbox_inches = "tight")

 指数積分 scipy.special.expi(x)

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

# 区間[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.savefig("expn.png", bbox_inches = "tight")

 指数積分 scipy.special.expinx)

 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\]
を計算します。

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.savefig("dawson.png", bbox_inches = "tight")

 Dawson integral (ドーソン積分)

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 のような特殊な記号を指定して広義積分を求めることもできます。