Scipyの数値積分のルーチンとしてGaussian quadratureを使うquadとかSimpson則simpsがあります。これの使い方自体はそこら中にドキュメントがあるかと思いますが、ここでちょっと定義が難しい関数を渡してみます。
言い訳免責事項
ここからは各種ドキュメントを斜め読みした浅い理解に基づいています。間違っている可能性も多々あります。
対象の関数
極限が収束する、あるいは解析的には問題なく計算できるが、素直にプログラムにするとエラーが出るような関数として、以下を例にしてみます。
f(x) = \frac{\sin(x)}{x}
こちらの関数について、$0$から$2\pi$までの積分
I = \int^{2\pi}_0\mathrm{d}x\frac{\sin(x)}{x}
を数値積分してみます。
問題点
これは解析的には $x\rightarrow 0$ で $f(x) \rightarrow 1$ になり、また積分もディレクレ積分としてあるように可能です。
ですがこれを何も考えず、
def f(x):
return np.sin(x)/x
と定義しますと、
>>> f(0.0)
RuntimeWarning: invalid value encountered in double_scalars
nan
となってしまいます。普通はかなり小さいxに対して1を返すような処理を差し込んだりするかと思いますが、今回はあえてこのまま突っ込んでみます。
答え合わせ
この関数の不定積分に対し、$\mathrm{Si}(z)$という特殊関数が定義されています。
http://mathworld.wolfram.com/SineIntegral.html
\mathrm{Si}(z) = \int^z_0\mathrm{d}t\frac{\sin(t)}{t}
これを用いると、うまく計算できると
I = \mathrm{Si}(2\pi) \approx 1.4181515761326284
となります。
ちなみにこの特殊関数はscipyでも提供されています。
from scipy.special import sici
from numpy import pi
si, ci = sici(2.0*pi)
print(si) # 1.4181515761326284
siciは$\mathrm{Si}(z)$のCosine版$\mathrm{Ci}(z) = \gamma + \log(x) + \int^z_0\frac{\cos(z) - 1}{z}$とのタプルとして返します。
環境
- Ubuntu 16.04
- python 3.5
- numpy 1.14.3
- scipy 1.0.0
Gaussian quadrature
何も考えずquadを使ってみます。
from scipy.integrate import quad
def f(x):
return np.sin(x)/x
I, _ = quad(f, 0.0, 2.0*pi)
print(I)
これを実行してみると・・・
1.4181515761326284
あれ、普通に行けた。しかも$ \mathrm{Si}(2\pi) \approx 1.4181515761326284 $ が表示した桁全部で一致しています。
もしかしたら $x = 0$ が積分域の端にあると積分点にならないのかもしれない。
積分区間を変えて
定義から分かるように、$f(x)$は偶関数なので、$[-2\pi, 2\pi]$の積分は$[0, 2\pi]$の積分の２倍です。ということで、
I, _ = quad(f, -2.0*pi, 2.0*pi)
print(I*0.5)
としてみると、
RuntimeWarning: invalid value encountered in double_scalars
return np.sin(x)/x
/usr/local/lib/python3.5/dist-packages/scipy/integrate/quadpack.py:364: IntegrationWarning: The occurrence of roundoff error is detected, which prevents
the requested tolerance from being achieved. The error may be
underestimated.
warnings.warn(msg, IntegrationWarning)
nan
となりました。
特異点ありGaussian quadrature
quad関数にはpointsというオプション引数があります。
scipyのquadのマニュアルから引用しますと：
points : (sequence of floats,ints), optional
A sequence of break points in the bounded integration interval where local difficulties of the integrand may occur (e.g., singularities, discontinuities). The sequence does not have to be sorted.
とあります。こいつに $x = 0$を渡してやってみようと思います。
まずは $[0, 2\pi]$の積分について。
I, _ = quad(f, 0.0, 2.0*pi, points=(0.0, ))
print(I)
/usr/local/lib/python3.5/dist-packages/scipy/integrate/quadpack.py:364: IntegrationWarning: Extremely bad integrand behavior occurs at some points of the
integration interval.
warnings.warn(msg, IntegrationWarning)
nan
一方で
I, _ = quad(f, -2.0*pi, 2.0*pi, points=(0.0, ))
print(I*0.5)
とすると何のエラーや警告なく
1.4181515761326284
といけます。
(追記：scipyを1.1.0にバージョンを上げると両方のケースで問題なく積分できます)
チェック
quadに渡したpointsは、そこで積分区間を分割するということらしいです。ということでどう分割されたかをチェックしてみます。
I, _, info, _ = quad(f, 0.0, 2.0*pi, points=(0.0, ), full_output=1)
print(info["pts"])
I, _, info = quad(f, -2.0*pi, 2.0*pi, points=(0.0, ), full_output=1)
print(info["pts"])
full_outputに0でない値を渡してやるとinfoに当たる値が帰ってきて、これで積分に関するいろいろな情報をチェックできます。ちなみに１回目と２回目で返り値を受け取る変数の数が違うのは、１回目の方はエラー情報を返してくるからです。
これを実行すると、
[0. 0. 6.28318531]
[-6.28318531 0. 6.28318531]
つまり、$[0, 2\pi]$ の積分で $x = 0$をpointsに渡すとエラーになるのは、 $[0, 0]$の積分が実行され、f(0.0)が評価されてしまったからかと思います。
Simpson則
これはあらかじめy=f(x)の値を用意しておく必要があります。
xs = np.arange(0.0, 2.0*pi, 0.0001)
ys = f(xs)
こうすると、ys[0]がnanになって、どうせ積分値もnanになります。そこで、$[-2\pi, 2\pi]$の区間の積分として、またdxをうまく調整して $x = 0$ を避けます。
xs = np.arange(-2.0*pi, 2.0*pi, 0.0001)
ys = f(xs)
print(len(xs))
I = simps(ys, xs)
print(I*0.5)
これを実行すると
125664
1.4181515763310302
125664個のサンプル数を持ってしてもquadほど正確ではないようです。(参考：$ \mathrm{Si}(2\pi) \approx 1.4181515761326284 $ )
ただsimpsはfunctionとして定義できないような測定データなどを渡せるという利点があるかと思います。
追記
この記事を書いたあと、ふとscipyを1.1.0にバージョンアップしてみたのですが、その結果
I, _, info = quad(f, 0.0, 2.0*pi, points=(0.0, ), full_output=1)
print(I)
print(info["pts"])
に対して
1.4181515761326284
[0. 6.28318531]
と、$[0,0]$となる積分区間がうまくスキップされるようになっていて、問題なく積分できるようになっていました。
参考