Scipyの数値積分の使い方

Scipyの数値積分のルーチンとしてGaussian quadratureを使う`quad`とかSimpson則`simps`があります。これの使い方自体はそこら中にドキュメントがあるかと思いますが、ここでちょっと定義が難しい関数を渡してみます。

# <strike>言い訳</strike>免責事項

ここからは各種ドキュメントを斜め読みした浅い理解に基づいています。間違っている可能性も多々あります。

# 対象の関数

極限が収束する、あるいは解析的には問題なく計算できるが、素直にプログラムにするとエラーが出るような関数として、以下を例にしてみます。

```math
f(x) = \frac{\sin(x)}{x}
```

![sinex_x.png](https://qiita-image-store.s3.amazonaws.com/0/258019/9efbdde1-f2bc-3fa2-f53e-718024e574f4.png)

こちらの関数について、$0$から$2\pi$までの積分

```math
I = \int^{2\pi}_0\mathrm{d}x\frac{\sin(x)}{x}
```
を数値積分してみます。

# 問題点

これは解析的には $x\rightarrow 0$ で $f(x) \rightarrow 1$ になり、また積分もディレクレ積分としてあるように可能です。
ですがこれを何も考えず、

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

と定義しますと、

```python
>>> f(0.0)
RuntimeWarning: invalid value encountered in double_scalars
nan
```

となってしまいます。普通はかなり小さい`x`に対して`1`を返すような処理を差し込んだりするかと思いますが、今回はあえてこのまま突っ込んでみます。

# 答え合わせ

この関数の不定積分に対し、$\mathrm{Si}(z)$という特殊関数が定義されています。
http://mathworld.wolfram.com/SineIntegral.html

```math
\mathrm{Si}(z) = \int^z_0\mathrm{d}t\frac{\sin(t)}{t}
```

これを用いると、うまく計算できると

```math
I = \mathrm{Si}(2\pi) \approx 1.4181515761326284
```

となります。

ちなみにこの特殊関数はscipyでも提供されています。

```python
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`を使ってみます。

```python
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]$の積分の２倍です。ということで、

```python
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]$の積分について。

```python
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
```
一方で

```python
I, _ = quad(f, -2.0*pi, 2.0*pi, points=(0.0, ))
print(I*0.5)
```
とすると何のエラーや警告なく

```
1.4181515761326284
```
といけます。

(追記：scipyを1.1.0にバージョンを上げると両方のケースで問題なく積分できます)



## チェック

`quad`に渡した`points`は、そこで積分区間を分割するということらしいです。ということでどう分割されたかをチェックしてみます。

```python
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)`の値を用意しておく必要があります。

```python
    xs = np.arange(0.0, 2.0*pi, 0.0001)
    ys = f(xs)
```
こうすると、`ys[0]`が`nan`になって、どうせ積分値も`nan`になります。そこで、$[-2\pi, 2\pi]$の区間の積分として、また`dx`をうまく調整して $x = 0$ を避けます。

```python
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にバージョンアップしてみたのですが、その結果

```python
    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]$となる積分区間がうまくスキップされるようになっていて、問題なく積分できるようになっていました。


# 参考

https://github.com/scipy/scipy/issues/7555