積分をPythonで理解する


積分もやはり、Sympy というライブラリが非常に便利です。

# Sympy を使った方法


```python
import sympy as sym
from sympy.plotting import plot
sym.init_printing(use_unicode=True)
%matplotlib inline
```

指定した文字をシンボル（変数を表す文字）として扱います。


```python
a, b, c, x, y = sym.symbols("a b c x y")
```

主要な数学関数を sympy からインポートしておきます。


```python
from sympy import sin, cos, tan, log, exp
```

## 積分


```python
# 積分の式
expr  = x ** a
integ = sym.Integral(expr, x)
print(integ)
integ
```

    Integral(x**a, x)






$$\int x^{a}\, dx$$




```python
# 積分の実行
integ.doit()
```




$$\begin{cases} \frac{x^{a + 1}}{a + 1} & \text{for}\: a \neq -1 \\\log{\left (x \right )} & \text{otherwise} \end{cases}$$




```python
# このように書いても同じ
sym.integrate(expr, x)
```




$$\begin{cases} \frac{x^{a + 1}}{a + 1} & \text{for}\: a \neq -1 \\\log{\left (x \right )} & \text{otherwise} \end{cases}$$




```python
# sym.Eq を使って等式として表記
eq = sym.Eq(sym.Integral(expr, x), sym.integrate(expr, x))
print(eq)
eq
```

    Eq(Integral(x**a, x), Piecewise((x**(a + 1)/(a + 1), Ne(a, -1)), (log(x), True)))





$$\int x^{a}\, dx = \begin{cases} \frac{x^{a + 1}}{a + 1} & \text{for}\: a \neq -1 \\\log{\left (x \right )} & \text{otherwise} \end{cases}$$



## いろんな積分の公式


```python
expr = 1/x
eq = sym.Eq(sym.Integral(expr, x), sym.integrate(expr, x))
print(eq)
eq
```

    Eq(Integral(1/x, x), log(x))





$$\int \frac{1}{x}\, dx = \log{\left (x \right )}$$




```python
expr = a ** x
eq = sym.Eq(sym.Integral(expr, x), sym.integrate(expr, x))
print(eq)
eq
```

    Eq(Integral(a**x, x), Piecewise((a**x/log(a), Ne(log(a), 0)), (x, True)))





$$\int a^{x}\, dx = \begin{cases} \frac{a^{x}}{\log{\left (a \right )}} & \text{for}\: \log{\left (a \right )} \neq 0 \\x & \text{otherwise} \end{cases}$$




```python
expr  = exp(x)
eq = sym.Eq(sym.Integral(expr, x), sym.integrate(expr, x))
print(eq)
eq
```

    Eq(Integral(exp(x), x), exp(x))





$$\int e^{x}\, dx = e^{x}$$




```python
expr = sin(x)
eq = sym.Eq(sym.Integral(expr, x), sym.integrate(expr, x))
print(eq)
eq
```

    Eq(Integral(sin(x), x), -cos(x))





$$\int \sin{\left (x \right )}\, dx = - \cos{\left (x \right )}$$




```python
expr = cos(x)
eq = sym.Eq(sym.Integral(expr, x), sym.integrate(expr, x))
print(eq)
eq
```

    Eq(Integral(cos(x), x), sin(x))





$$\int \cos{\left (x \right )}\, dx = \sin{\left (x \right )}$$




```python
expr = tan(x)
eq = sym.Eq(sym.Integral(expr, x), sym.integrate(expr, x))
print(eq)
eq
```

    Eq(Integral(tan(x), x), -log(cos(x)))





$$\int \tan{\left (x \right )}\, dx = - \log{\left (\cos{\left (x \right )} \right )}$$




```python
expr = log(x)
eq = sym.Eq(sym.Integral(expr, x), sym.integrate(expr, x))
print(eq)
eq
```

    Eq(Integral(log(x), x), x*log(x) - x)





$$\int \log{\left (x \right )}\, dx = x \log{\left (x \right )} - x$$




```python
expr = 1/cos(x)**2
eq = sym.Eq(sym.Integral(expr, x), sym.integrate(expr, x))
print(eq)
eq
```

    Eq(Integral(cos(x)**(-2), x), sin(x)/cos(x))





$$\int \frac{1}{\cos^{2}{\left (x \right )}}\, dx = \frac{\sin{\left (x \right )}}{\cos{\left (x \right )}}$$




```python
expr = 1/sin(x)**2
eq = sym.Eq(sym.Integral(expr, x), sym.integrate(expr, x))
print(eq)
eq
```

    Eq(Integral(sin(x)**(-2), x), -cos(x)/sin(x))





$$\int \frac{1}{\sin^{2}{\left (x \right )}}\, dx = - \frac{\cos{\left (x \right )}}{\sin{\left (x \right )}}$$




```python
expr = (x-a)**b
eq = sym.Eq(sym.Integral(expr, x), sym.integrate(expr, x))
print(eq)
eq
```

    Eq(Integral((-a + x)**b, x), Piecewise(((-a + x)**(b + 1)/(b + 1), Ne(b, -1)), (log(-a + x), True)))





$$\int \left(- a + x\right)^{b}\, dx = \begin{cases} \frac{\left(- a + x\right)^{b + 1}}{b + 1} & \text{for}\: b \neq -1 \\\log{\left (- a + x \right )} & \text{otherwise} \end{cases}$$




```python
expr = 1/sin(x)
eq = sym.Eq(sym.Integral(expr, x), sym.integrate(expr, x))
print(eq)
eq
```

    Eq(Integral(1/sin(x), x), log(cos(x) - 1)/2 - log(cos(x) + 1)/2)





$$\int \frac{1}{\sin{\left (x \right )}}\, dx = \frac{\log{\left (\cos{\left (x \right )} - 1 \right )}}{2} - \frac{\log{\left (\cos{\left (x \right )} + 1 \right )}}{2}$$



# 定積分


```python
expr = (x-a)*(b-x)
eq = sym.Eq(sym.Integral(expr, (x, a, b)), sym.integrate(expr, (x, a, b))).factor()
print(eq)
eq
```

    Eq(-Integral((-a + x)*(-b + x), (x, a, b)), -(a - b)**3/6)





$$- \int_{a}^{b} \left(- a + x\right) \left(- b + x\right)\, dx = - \frac{\left(a - b\right)^{3}}{6}$$




```python
expr = x/(x**2 + 1)
eq = sym.Eq(sym.Integral(expr, (x, 1, 2)), sym.integrate(expr, (x, 1, 2)))
print(eq)
eq
```

    Eq(Integral(x/(x**2 + 1), (x, 1, 2)), -log(2)/2 + log(5)/2)





$$\int_{1}^{2} \frac{x}{x^{2} + 1}\, dx = - \frac{\log{\left (2 \right )}}{2} + \frac{\log{\left (5 \right )}}{2}$$



Sympy ではこれ以上展開してくれる（小数表示にしてくれる）方法が見つからなかったので、とりあえず numpyなどから関数を借りてきて次のように計算してみました。


```python
from numpy import log
-log(2)/2 + log(5)/2
```




$$0.4581453659370775$$



## 円周率

円周率は次のような積分で表されます。


```python
expr = 4/(x**2 + 1)
eq = sym.Eq(sym.Integral(expr, (x, 0, 1)), sym.integrate(expr, (x, 0, 1)))
print(eq)
eq
```

    Eq(Integral(4/(x**2 + 1), (x, 0, 1)), pi)





$$\int_{0}^{1} \frac{4}{x^{2} + 1}\, dx = \pi$$



# Scipy.integrate を使った方法


```python
from scipy import integrate
```

積分したい式を関数で用意します。


```python
func = lambda x: 4 / (x ** 2 + 1)
```

integrate.quad で、積分と、推定誤差が得られます。


```python
integrate.quad(func, 0, 1)
```




$$\left ( 3.1415926535897936, \quad 3.4878684980086326e-14\right )$$



## 台形則とシンプソン則

台形則やシンプソン則により積分するには、あらかじめ刻み幅を指定しておく必要があります。


```python
import numpy as np
x=np.linspace(0,1,5)
```


```python
integrate.trapz(func(x), x) # 台形則
```




$$3.131176470588236$$




```python
integrate.simps(func(x), x) # シンプソン則
```




$$3.1415686274509804$$


# 課題

* 積分を数値的に求める区分求積法について、長方形近似、台形則、シンプソン則の違いを説明してください。

* 上記の Integral(x/(x**2 + 1), (x, 1, 2)) の式が解析的に解けるのは上に示した通りですが、これを数値的に解く以下の３つのプログラムを Python で作成してください。ただし、Sympy も Scipy.integrate も使わないこと。インターネットで検索して見つけたプログラムを適宜改変して使っても良いですが、その場合は出典のURLを明記してください。
  * 長方形近似で解くプログラム
  * 台形則で解くプログラム
  * シンプソン則で解くプログラム