数値積分¶
SciPy を利用して、何か適当な実数値関数の定積分を計算する手順をいくつか記す。まずは区分求積法の計算法の代表例である台形公式と Simpson の公式を実装したと思われるSciPy の関数を紹介し、次に使い勝手のよい定積分計算関数を紹介する。
台形公式¶
関数 scipy.integrate.trapz を利用して、実数値関数の台形公式による定積分の近似値を得ることができる。手順は次のようなものだ。
被積分関数
fを Python の関数として定義または用意する。積分区間上の列
xを用意する。列
xに対し、関数を評価した値の列yを作る。scipy.integrate.trapzを呼び出す。
以下にコード例を示す。
#!/usr/bin/env python
"""integrate.py: Demonstrate some SciPy integration functions.
"""
from scipy.integrate import trapz, simps, quad
import numpy as np
# pylint: disable=invalid-name
def f(t):
"""Integrand function."""
return np.sin(t)
# Define sample points.
x = np.linspace(0, np.pi, num=100, endpoint=True)
y = f(x)
# Compute the integral by using trapezoid rule.
print("Trapezoid method: ", trapz(y, x))
\(\displaystyle \int_{0}^{\pi} \sin x\,\dd{x}\) という簡単な例で面白くないが、見ていこう。
まずは NumPy の関数 linspace を用いて、積分区間 \({[0, \pi]}\) を百等分した列 x を用意する。これはコード的には np.ndarray 型なので、単に {y =
np.sin(x)} とすることで関数値の列 y が得られる。これで計算の準備が整った。直ちに両者を関数 trapz の引数に与えて呼び出す。
実行結果は関連定積分関数の結果とまとめて後述する。
Simpson の公式¶
関数 scipy.integrate.simps を利用すると、前項同様の方法で定積分を計算できる。コードを、台形公式との違いのみを示す。
# Compute the integral by using Simpsons's rule.
print("Simpson's method: ", simps(y, x))
実行結果は関連定積分関数の結果とまとめて後述する。
汎用目的積分関数¶
関数 scipy.integrate.quad を利用すると、今度はサンプリングの手間を省いて定積分を計算することになる。これは、被積分関数本体と、積分区間から定積分を計算する。
以下にコード例を示す。
#!/usr/bin/env python
"""integrate.py: Demonstrate some SciPy integration functions.
"""
from scipy.integrate import trapz, simps, quad
import numpy as np
# pylint: disable=invalid-name
def f(t):
"""Integrand function."""
return np.sin(t)
# Define sample points.
x = np.linspace(0, np.pi, num=100, endpoint=True)
y = f(x)
# Compute the integral by using QUADPACK (Fortran).
# interval of integration; 0 to pi.
I3 = quad(f, x[0], x[-1])
print("QUADPACK: ", I3)
実行結果は関連定積分関数の結果とまとめて後述する。
この関数はどういう求積法を実装しているのか。ドキュメントによると <a technique from the Fortran library QUADPACK> を利用したとあるので、かなり専門的なのかもしれない。当関数には実際は大量の引数が存在しており、求積法に対して微細な指示ができそうだ。例えば今回の例では有限区間の定積分計算を試したが、値が求まる限り無限区間も計算可能のようだ。
実行結果¶
前述までの計算結果を一挙に記す。すべて \({\displaystyle \int_{0}^{\pi} \sin x\,\dx}\) である。
Trapezoid method: 1.99983216389
Simpson's method: 1.99999996902
QUADPACK: (2.0, 2.220446049250313e-14)
この例では関数 quad の精度が圧倒的によいようだ。
