なんかscipyにあるquadを使うと、Gauss求積法による数値積分ができるそうで、しかも無限区間も可能なようなので、試してみます。

環境

scipy 1.1.0
python 3.5.2

あと、基本的に以下のモジュールをimportしています。

from scipy.integrate import quad
import numpy as np

ちょろい例

有理関数

\int_{-\infty}^\infty \frac{\mathrm{d}x}{\left(x^2 + 1\right)^{n + 1}} = \frac{\pi(2n)!}{2^{2n}(n!)^2}

これの $n=3$ のときをscipyで計算してみます。

    n = 3
    def f(x):
        return 1.0/(x**2 + 1)**(n+1)
    I, _ = quad(f, -np.inf, np.inf)

    factorial = lambda n: np.prod(np.arange(1, n + 1))
    ref = np.pi*factorial(2*n)/(2**(2*n))/((factorial(n))**2)

    print("reference:", ref)
    print("QUAD:", I)
    print("Error:", np.abs(I - ref))

Iが数値積分の結果、refが解析解の参考値です。

これを実行してみると・・・

reference: 0.9817477042468103
QUAD: 0.9817477042468103
Error: 0.0

完璧ですね。

対数

\int^\infty_0\mathrm{d}x\frac{\log(x)}{(x + a)^2} = \frac{\log(a)}{a}

被積分関数は$ x = 0 $ で発散します。簡単のため、 $a=e$ (自然対数の底)としてみます。

    a = np.e
    def f(x):
        return np.log(x)/((x + a)**2)
    ref = np.log(a)/a
    I, _ = quad(f, 0, np.inf)
    print("reference:", ref)
    print("QUAD:", I)
    print("Error:", np.abs(I - ref))

これを実行すると

reference: 0.36787944117144233
QUAD: 0.36787944117149923
Error: 5.689893001203927e-14

少し誤差が出てきましたが、倍精度で計算していることを考えると十分な精度でしょう。

振動Gaussian

三角関数とGaussianの積です。なんか正式な名前があった気がします。

\int^\infty_{-\infty}\mathrm{d}x \exp( - x^2)\cos(2bx) = \sqrt{\pi}\exp(-b^2)

    b = 1.0
    def f(x): 
        return np.exp( -x*x )*np.cos(2.0*b*x) 
    ref = np.sqrt(np.pi)*np.exp(-b**2)
    I, _ = quad(f, -np.inf, np.inf)
    print("reference:", ref)
    print("QUAD:", I)
    print("Error:", np.abs(I - ref))

これを実行すると

reference: 0.6520493321732922
QUAD: 0.6520493321732844
Error: 7.771561172376096e-15

これも十分な精度がでているようです。

手こずった例

有理関数その２

Cauchyの主値積分の例題です。
http://eman-physics.net/math/imaginary12.html

PV\int^\infty_{^-\infty}\frac{\mathrm{d}x}{x^3 - x^2 + x - 1} = -\frac{\pi}{2}

この被積分関数は$x = 1$で発散します。複素関数論では、この点の周りを微小円で回避して、その微小円の半径を $+0$ にする極限を取るというストーリーになるかと思います。数値積分のquadでもそれに倣ってみます。
関数と参考値。

    def f(x):
        return 1.0/(x**3 - x**2 + x - 1)
    ref = - np.pi/2.0

これを素直に、 $x = 1$を回避しつつquadを使ってみます。

    I, _ = quad(f, -np.inf, np.inf, points=[1.0, ])

これを実行すると

ValueError: Infinity inputs cannot be used with break points.

残念。ということで手で分割してみます。

    I1, _ = quad(f, -np.inf, 1.0)
    I2, _ = quad(f, 1.0, np.inf)
    I = I1 + I2

    print("reference:", ref)
    print("QUAD:", I)
    print("Error:", np.abs(I - ref))

これを実行すると

/usr/local/lib/python3.5/dist-packages/scipy/integrate/quadpack.py:385: IntegrationWarning: Extremely bad integrand behavior occurs at some points of the
  integration interval.
  warnings.warn(msg, IntegrationWarning)
reference: -1.5707963267948966
QUAD: -2.164415573210107
Error: 0.5936192464152104

う〜ん。では微小円に相当することをやってみます。

    eps = 1e-5
    I1, _ = quad(f, -np.inf, 1.0 - eps)
    I2, _ = quad(f, 1.0 + eps, np.inf)
    I = I1 + I2

reference: -1.5707963267948966
QUAD: -1.5707863268002527
Error: 9.999994643905552e-06

epsのオーダーを考えると、良い精度ではないかと思います。
ちなみにいろいろやってみましたが、他のオプションを変えない場合 eps = 1e-10 あたりでErrorが一番小さく、1.807e-07 程度になりました。

ディリクレ積分

前回のScipyの数値積分の使い方の例では、区間 $[0, 2\pi]$ なのでかなり精度が高かったようです。今回は $[0, \infty]$ でやってみます。

\int^\infty_0\mathrm{d}x\frac{\sin(x)}{x} = \lim_{x\rightarrow\infty}\mathrm{Si}(x) = \frac{\pi}{2}

いろいろ試行錯誤した結果をそのまま貼り付けます。けっしてめんどくさくなったとかではないです。
full_outputオプションが使われているのは、いろいろ試していた痕跡で、特に意味はないです。

from scipy.integrate import quad
from scipy.special import sici
import numpy as np

def f(x):
    return np.sin(x)/x

def main():
    si, ci = sici(np.inf)
    print(si)
    print("-----------")
    I, _, info, msg = quad(f, 0.0, np.inf, full_output=1)
    print(I)
    print("-----------")
    I, _, info, msg = quad(f, 0.0, np.inf, limit=1000, full_output=1)
    print(I)
    print("-----------")
    I, _ = quad(f, 0.0, 1000000)
    print(I)
    print("-----------")
    I, _ = quad(f, 0.0, 1000000, limit=10000)
    print(I)
    print("-----------")
    I, _ = quad(f, 0.0, 1000000, limit=300000)
    print(I)
    print(np.abs(si - I))

if __name__ == '__main__':
    main()

これを実行すると

1.5707963267948966
-----------
2.2478679634689205
-----------
1.840170585978155
-----------
/usr/local/lib/python3.5/dist-packages/scipy/integrate/quadpack.py:385: IntegrationWarning: The maximum number of subdivisions (50) has been achieved.
  If increasing the limit yields no improvement it is advised to analyze 
  the integrand in order to determine the difficulties.  If the position of a 
  local difficulty can be determined (singularity, discontinuity) one will 
  probably gain from splitting up the interval and calling the integrator 
  on the subranges.  Perhaps a special-purpose integrator should be used.
  warnings.warn(msg, IntegrationWarning)
-3.150746880828421
-----------
/usr/local/lib/python3.5/dist-packages/scipy/integrate/quadpack.py:385: IntegrationWarning: The integral is probably divergent, or slowly convergent.
  warnings.warn(msg, IntegrationWarning)
0.2795622380611621
-----------
1.5707953900431475
9.367517490588284e-07

簡単にまとめると、

普通に区間 $[0, \infty]$で quadを使うと0.7程度ずれてしまいます。
区間 $[0, \infty]$でlimit（区間分割の上限値）を増やすと若干改善しますが、それでも0.3程度ずれます。
区間を $[0, 1000000]$とすると、limitが足りず精度が出てないという警告がでます。結果は大きくずれています。
そこでlimit=10000とすると、収束が遅いという警告がでます。結果の誤差は少しはマシですが、かなり大きいです。
limit=300000にするとようやく誤差が1e-6程度になりました。

なおlimitを大きくすると時間もかかります。

$ time python dirichlet_inf.py
（標準出力、中略）
python dirichlet_inf.py  11.34s user 0.14s system 100% cpu 11.454 total

これまでほぼ一瞬で終わっていたquadですが、limit=10000で1秒弱、limit=300000とすると10秒程度かかるようになりました。

今回は絶対収束しない積分で収束が遅く、このような工夫が必要となったようです。

参考

積分例題はここを参考にしました。
http://www.th.phys.titech.ac.jp/~muto/lectures/Amath06/ex_chap12.pdf

ディリクレ積分が難しいという話
https://stackoverflow.com/questions/46509934/calculating-dirichlet-integral-in-python

感想

quadに渡したfは被積分関数の表面的な定義をそのまま渡しているだけで、$0$、特異点や $\infty$などでの挙動を一切教えていないのに、これほどの精度が出るとは思っていませんでした。
手こずった例も、結局関数そのものを工夫することなく高精度を目指せるようです。

ふと、複素関数論の授業のレポート問題で、めんどくさくなってMaximaに投げたものの、途中経過が一切ないので泣く泣く自力で頑張ったことを思い出しました。

その他

僕はGaussian quadratureリテラシーは高くないので、weightの使い方はわかりません。

scipyのquadベンチマーク

環境

ちょろい例

有理関数

対数

振動Gaussian

手こずった例

有理関数その２

ディリクレ積分

参考

感想

その他