数値積分

当サイトではアフィリエイトプログラムを利用して商品を紹介しています。

この記事では、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$ です)。
 
方法の1つは 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

もう1つの方法は 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 台形公式 (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) にしたがって数値積分を実行します。

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 フレネル積分のグラフ

指数積分

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.expi(x)
 
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.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\]
を計算します。

# 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()

Dawson integral (ドーソン積分)

【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)
\[x\log x-x\]

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

# In[2]

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

display(integ)
\[-1+2\log(2)\]

そのままでは、$\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)
\[\pi\]

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

# In[5]

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

display(integ)
\[\log{\left (\cos(x) \right )} + \frac{1}{2 \cos^{2}{\left (x \right )}}\]

コメント

  1. HNaito より:

    下記は誤植と思われますので、ご確認ください。
     
    DOUBLE_INTEGRAL In[1] プログラムのコメントで、# yを0からpiまで → fを0から1まで
     
    DOUBLE_INTEGRAL In[1] と In[2] で、積分範囲をlambda関数で指定していますが、そのまま 0 と 1 を書いても同じ結果が得られました。
     
    EI_INTEGRAL In[2] プログラムの下の文で、scipy.special.exp1(x, n) → scipy.special.exp1(x)
    SYMPY_INTEGRAL In[1] プログラムの上の文で、∫ log x → ∫ log x dx
    SYMPY_INTEGRAL In[3] プログラムの下の文で、∫1/(1+x^2) → ∫ 1/(1+x^2) dx
    SYMPY_INTEGRAL In[5] プログラムの実行結果が、log(cos(x)) + (1/(2cos^2(x)) となりました。

  2. あとりえこばと より:

    【AI連載小説】科学とコードの交差点(3)
     
    開誠は微分の学習に続き、積分についても理解を深めるために講義に参加していた。教室の雰囲気は静かで、熱心な学生たちが座席に座り、期待と興奮を胸に授業の始まりを待っていた。教室の前に立つ教授は、物理学における積分の重要性について導入を行っていた。
    「皆さん、おはようございます。今日からは微分に続き、積分について学んでいきます。微分と積分は対を成すもので、物理学においては不可欠なツールです」
    開誠は興奮を抑えきれない様子でメモを取りながら、教授の話に耳を傾けていた。
    「微分が関数の傾きや変化率を表すのに対し、積分は関数の面積や量を表します。これらは物理学において、例えば加速度から速度、力から仕事といった概念を導き出すのに使われます」
    白板には積分の記号や基本的な公式が書かれ、教授は具体的な例を挙げながら解説していった。
    「具体的な積分の計算にはいくつかの手法があります。不定積分、定積分、そして数値積分などです」
    授業は進むにつれて、開誠と同様に熱心な学生たちは積分の理論に没頭し、質問や議論が交わされていた。開誠も積分の概念がより深く理解できるように、積極的に参加し、教授の指導のもと新しい知識を吸収していった。

    積分の計算をより効率的に行うため、開誠はPythonで使える積分計算ライブラリを見つけるためにオンラインで検索を始めた。彼は積分に関する様々なライブラリやツールを調査し、どれが彼の目的に最適かを見極めようとしていた。ブラウザの検索バーに「Python integration library」を入力し、検索結果を一つ一つ見ていった。

    検索結果:
    1. “SciPy” – 数値計算、積分計算もサポート
    2. “SymPy” – 記号計算ライブラリ、不定積分なども扱える
    3. “NumPy” – 数値計算の基本ライブラリ、積分関数も提供
    4. “Quadpy” – 数値積分のためのライブラリ、多様な積分手法が利用可能
    5. “GSL” – GNU Scientific Library、高度な数値計算や積分が可能
    6. “Scikit-learn” – 機械学習ライブラリ、積分も一部サポート
     
    開誠はそれぞれのライブラリが提供する機能や利点を比較し、自分の問題に最も適したものを見つけようとしていた。特に数値計算ライブラリであるSciPyやQuadpyは、彼が物理学の問題に取り組む上で有望な選択肢となりそうだった。