Pythonで積分

※ この記事はnoteをやっていた頃に投稿した記事です。こちらのブログに移しました。
（note投稿日：2022年9月9日）

初めに
integrate.quadについて
1. 引数について
2. 返り値について
実際に使ってみる
まとめ

初めに

SciPy(Pythonのライブラリ)を使って積分の計算ができるらしい

ドキュメントはこちら
「SciPy」から「integrate.quad」を利用して積分を計算する

integrate.quadについて

引数について

quad(func, a, b, args=( ), full_output=0, epsabs=1.49e-08, epsrel=1.49e-08, limit=50, points=None, weight=None, wvar=None, wopts=None, maxp1=50, limlst=50)
ドキュメントまたはソースコードより

引数	説明
func	関数
a	積分区間の下端
b	積分区間の上端
args	パラメータ

※ただし、argsは積分したい関数にパラメータを持たせたいときに使う

返り値について

(y, abserr)の形で返ってくる

返り値	型	説明
y	float	積分の値
abserr	float	絶対誤差の推定値

とりあえずこれだけ知っていれば計算できそう。後は必要になったら調べるとしよう。

実際に使ってみる

簡単な例

ドキュメントと同じ例だが、

$$
\int_{0}^4x^2dx
$$
を計算してみる。

コード

from scipy import integrate

y = lambda x: x**2

result, absolute_error = integrate.quad(y,0,4)

print(result)
print(absolute_error)

出力

21.333333333333332 //result
2.3684757858670003e-13 //absolute_error

実際

$$
\int_{0}^4x^2dx = \frac{64}{3} = 21.333333333333…
$$

パラメータを加えてみた例

次にパラメータを加えてみる。

$$
\int_{0}^{4}ax^2dx
$$

コード

from scipy import integrate

y = lambda x: x**2

result, absolute_error = integrate.quad(y,0,4)

print(result)
print(absolute_error)```python
from scipy import integrate

y = lambda x: x**2

result, absolute_error = integrate.quad(y,0,4)

print(result)
print(absolute_error)
```from scipy import integrate

y = lambda x: x**2

result, absolute_error = integrate.quad(y,0,4)

print(result)
print(absolute_error)

出力

21.333333333333332 //result
2.3684757858670003e-13 //absolute_error出力

```
21.333333333333332 //result
2.3684757858670003e-13 //absolute_error
```

広義積分1

$$
\int_{-\infty}^{\infty}\frac{1}{\sqrt{2\pi}}\exp\left(-\frac{x^2}{2}\right)
$$

を計算してみる。(実際の積分結果は1になる)

from scipy import integrate
import numpy as np

y = lambda x: (1/(np.sqrt(2*np.pi)))*np.exp(-(x**2)/2)

result, absolute_error = integrate.quad(y, -np.inf, np.inf)

print(result)
print(absolute_error)

出力

0.9999999999999997 //result
1.0178191380347127e-08 //absolute_error

広義積分2

$$
\int_{0}^{\infty}xdx
$$

この場合はどのような結果になるだろう？

from scipy import integrate
import numpy as np

y = lambda x: x

result, absolute_error = integrate.quad(y,0,np.inf)

print(result)
print(absolute_error)

出力

0.4999999961769933 //result
5.7336234760563265e-06 //absolute_error

???

あ、よく読んでみるとこんなのが書いてあった！

IntegrationWarning: The integral is probably divergent, or slowly convergent.
  result, absolute_error = integrate.quad(y,0,np.inf)

つまり、「積分が発散しているかゆっくり収束しているかと思われる」的なことを言われている。

広義積分3

$$
\int_{0}^{1}\frac{1}{x}dx
$$

この場合はどうなるだろう。

from scipy import integrate
import numpy as np

y = lambda x: 1/x

result, absolute_error = integrate.quad(y,0,1)

print(result)
print(absolute_error)

出力

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.
  result, absolute_error = integrate.quad(y,0,1)
41.67684067538801 //result
9.35056037314051 //absolute_error

Warningなげー、、、

とりあえず、「細分の最大の数50が達成された。limitをあげても改善されなければ、、、」的な感じで書いているのでlimitを上げてみる。

(デフォルトは50みたい)

from scipy import integrate

y = lambda x: 1/x

result, absolute_error = integrate.quad(y,0,1,limit=10000)

print(result)
print(absolute_error)

limit=10000まで上げると

IntegrationWarning: Extremely bad integrand behavior occurs at some points of the
  integration interval.
  result, absolute_error = integrate.quad(y,0,1,limit=10000)
709.8707227351608
9.35056037314736

Warningが変わった。

積分区間のいくつかの点で良くない振る舞いが起きている、的な。

実際、

$$
\int_{0}^{1}\frac{1}{x}dx = \log|x| + C
$$

$x=0$ では定義されないので、確かにそうだなて感じ。