【SciPy】数値積分
SciPy の積分パッケージ scipy.integrate.quad() は Fortran のライブラリ QUADPACK を使って、ガウスの数値積分 (Gaussian quadrature) を実行します。必須引数は被積分関数 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.0points には複数の値を渡せるので、積分範囲が複数の特異点を含む場合は、特異点のリストを渡してください。
重積分
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.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) にしたがって数値積分を実行します。
シンプソンの公式は曲線を 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 による不定積分を求めることができます。
定積分を計算する場合は次のように記述します。
例として、対数関数 $\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)