Pythonで数学の勉強：微分積分学

# 微分積分学

[数の扱い方](http://qiita.com/tibigame/items/e6f3e7a2033dc5bdec5d)
[数式と関数](http://qiita.com/tibigame/items/61cecf86fc978628bfee)
[グラフを描く](http://qiita.com/tibigame/items/0e48673179aeeab55a43) (本文中のplot関数はここで作ってます)
と進めてきてようやく数学的らしいことを議論できる土台が整ったので微積をやってみる。

### import

```py3
import sympy as sym
sym.init_printing()
Pi = sym.S.Pi # 円周率
E = sym.S.Exp1 # 自然対数の底
I = sym.S.ImaginaryUnit # 虚数単位
oo = sym.oo # 無限大

# 使用する変数の定義(小文字1文字は全てシンボルとする)
(a,b,c,d,e,f,g,h,i,j,k,l,m,n,o,p,q,r,s,t,u,v,w,x,y,z) = sym.symbols('a b c d e f g h i j k l m n o p q r s t u v w x y z')

import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline

from mpl_toolkits.mplot3d import Axes3D
```

## 極限

sympyでは無限大としてoo = sym.ooが使われる。

### 右側極限
$$\lim_{n\to +0}\frac{1}{n} = \infty$$

```py3
sym.limit(1/x, x, 0)
```

### 左側極限
$$\lim_{n\to -0}\frac{1}{n} = -\infty$$

```py3
sym.limit(1/x, x, 0, dir='-')
```

### 2変数関数

```py3
sym.limit((x**3+y**3)/(x-y), x, y)
```
$\infty$sign$(y^3)$
これはおそらく$x\to y+0$を計算している。

#### 問
次の極限を計算せよ。
$$\lim_{x\to \infty}\frac{2^{(3x-1)}}{(x^4e^{x})^2}$$

```py3
Y = x**4 * E**x
sym.limit(2**(3*x-1)/(Y**2), x, oo)
```
$$\infty$$

```py3
plot(2**(3*x-1)/(Y**2), (650, 730, 1), (0, 100), close=True)
```
![hiki.png](https://qiita-image-store.s3.amazonaws.com/0/133885/627335ec-cabc-1c18-d24c-ae9fb420631e.png)
本質的には$2^3>e^2$が効いている。

## テイラー展開

関数$f(x)$の$x=a$でのテイラー展開は
$$\sum_{n=0}^{\infty}\frac{f^{(n)}(a)}{n!}(x-a)^n$$

$f(x)=\frac{\sin(x)}{x}$をテイラー展開してみる。

```py3
sym.series(sym.sin(x)/x, x, x0=0, n=15)
```
$$1 - \\frac{x^{2}}{6} + \\frac{x^{4}}{120} - \\frac{x^{6}}{5040} + \\frac{x^{8}}{362880} - \\frac{x^{10}}{39916800} + \\frac{x^{12}}{6227020800} - \\frac{x^{14}}{1307674368000} + \\mathcal{O}\\left(x^{15}\\right)$$

```py3
plot(sym.sin(x)/x, (-10, 10, 0.05), (-0.5, 1), cn=10)
plot(1 - x**2/6 + x**4/120 - x**6/5040 + x**8/362880
     - x**10/39916800 + x**12/6227020800 - x**14/1307674368000, (-10, 10, 0.05), (-0.5, 1), close=True)
```
![hiki2.png](https://qiita-image-store.s3.amazonaws.com/0/133885/4e128b60-ad8f-4495-9d34-6d2e9bcb8494.png)
原点から離れるほどこの近似は役に立たない。

## 微分

微分とは微小な変分を調べるもので接戦の傾きなど関数の計量に役立つ。
sym.diff(関数, 微分するシンボル(いくつでも))で求められる。

```py3
sym.diff((sym.tan(x)*x)/(x+sym.log(x)), x)
```
$$\\frac{x \\left(-1 - \\frac{1}{x}\\right) \\tan{\\left (x \\right )}}{\\left(x + \\log{\\left (x \\right )}\\right)^{2}} + \\frac{x \\left(\\tan^{2}{\\left (x \\right )} + 1\\right)}{x + \\log{\\left (x \\right )}} + \\frac{\\tan{\\left (x \\right )}}{x + \\log{\\left (x \\right )}}$$
微分は積、分数、合成関数の微分なども公式で計算できるので、
微分できる関数の組み合わせだと微分可能だと思っていい。

```py3
sym.diff(x**4+sym.exp(2*x), x, 4) # 4階微分
```
$$8(2e^{2x}+3)$$

```py3
sym.diff(x**3 * y**5 * z, x, y, 2, z) # xで微分→yで2階微分→zで微分
```
$$60x^2y^3$$

### 偏微分の記号のまま残す

```py3
XX = sym.Derivative((x**5 * y**2 * z**4), x, y, z)
XX
```
$$\\frac{\\partial^{3}}{\\partial x\\partial y\\partial z} \\left(x^{5} y^{2} z^{4}\\right)$$

sympyでは大文字の関数は遅延評価用。(数式表示用)

### 遅延評価

```py3
XX.doit()
```
$$40x^4yz^3$$

## 積分

積分とは面積や体積を計量するもの。

### 不定積分

```py3
sym.integrate(y/x, x)
```

$$y\log(x)$$


### 定積分

$$\int_1^{\infty}\frac{1}{x^2}dx$$

```py3
sym.integrate(1/x**2, (x, 1, oo))
```
1

積分は微分と違って求まるとは限らない。

```py3
sym.integrate(sym.sin(x)/(1+sym.exp(x)), (x, 0, Pi))
```

$$\\int_{0}^{\\pi} \\frac{\\sin{\\left (x \\right )}}{e^{x} + 1}\\, dx$$

そこで数値積分する。

```py3
X = lambda x: np.sin(x)/(1+np.exp(x))

from scipy import integrate
print(integrate.quad(X, 0, np.pi))
```
(0.38522681353410615, 1.0755139989210769e-14) # 2つ目は推定誤差

### サイクロイドの面積

```py3
X = (t-sym.sin(t))
Y = (1-sym.cos(t))
plotp([X, Y], (0, 2* np.pi, 400), [0, 2 * np.pi, 0, np.pi], figs=(6, 3), close=True)
```
![hiki3.png](https://qiita-image-store.s3.amazonaws.com/0/133885/a601a7cf-c689-14d9-aa0b-d09379aecaa3.png)

```py3
sym.integrate(Y * sym.diff(X, t), (t, 0, 2*Pi))
```
$3\pi$

## 特殊関数
### ディラックのデルタ関数

```py3
D = sym.DiracDelta(x)
D
```
$$\delta(x)$$

```py3
sym.integrate(D, (x, -1, 1))
```
$1$

$x=0$以外では0だが0を含む区間で積分すると1になる。
$x=0$の1点に面積1が集中しているイメージ。

### ヘヴィサイドのステップ関数

```py3
H = sym.Heaviside(x)
plot(H,(-1, 1, 0.01), (-0.5, 1.5), close=True)
```
![hiki4.png](https://qiita-image-store.s3.amazonaws.com/0/133885/36810cbb-c263-0b98-02a9-0352f06c948d.png)
デルタの積分がステップ、ステップの微分がデルタである。
$x<0$を考えたくないからステップかけるとか…。

## 重積分
$z=\sqrt{1-x^2-y^2}$と$z=0$で囲まれる部分の体積を求めてみよう。
言わずと知れた半球である。

```math
V=\int\!\!\!\int_D \!\sqrt{1-x^2-y^2} dxdy
```

$D$の第一象限は$y=\sqrt{1-x^2}$だから、

```math
V=4\int_0^1\!\!\!\int_0^{\sqrt{1-x^2}} \!\!\!\sqrt{1-x^2-y^2} dydx
```

この積分は普通に求まらないから数値積分してみる。

```py3
print(2*np.pi/3) # 理論値
integrate.dblquad(lambda x, y: 4*np.sqrt(1-x**2-y**2), 0, 1, lambda x: 0, lambda x: np.sqrt(1-x**2))
```
2.0943951023931953
(2.0943951023931957,3.5335157022586827e−10)

 さて$x=r\cos{\theta}, y=r\sin{\theta}$とおくと、

```math
V=\int_0^{2\pi} \!\!\! \int_0^1 \!\!\sqrt{1-r^2}\ \begin{vmatrix}
\frac{\partial x}{\partial r} & \frac{\partial x}{\partial \theta} \\
\frac{\partial y}{\partial r} & \frac{\partial y}{\partial \theta} \\
\end{vmatrix} drd\theta
```

```py3
X = r * sym.cos(t)
Y = r * sym.sin(t)
Z = sym.sqrt(1-x**2-y**2).subs([(x, X), (y, Y)]).trigsimp() # 簡略化は必要(1行ずつやった方がいい)
Jacobian = sym.Matrix([[sym.diff(X, r), sym.diff(X, t)], [sym.diff(Y, r), sym.diff(Y, t)]])
sym.integrate(sym.integrate(Z * Jacobian.det().trigsimp(), (r, 0, 1)), (t, 0, 2*Pi))
```

$$\frac{2\pi}{3}$$

見方を変えれば$V$の全領域を$1$で重積分しても同じだから
$$V=\int_V dxdydz$$
$x=r\sin{\theta}\cos{\varphi}, y=r\sin{\theta}\sin{\varphi}, z=r\cos{\theta}$とおくと、

```math
V=2\int_0^{\!\frac{\pi}{2}}\!\!\!\int_0^{\pi}\!\!\!\int_0^1
\begin{vmatrix}
\frac{\partial x}{\partial r} & \frac{\partial x}{\partial \theta} & \frac{\partial x}{\partial \varphi} \\
\frac{\partial y}{\partial r} & \frac{\partial y}{\partial \theta} & \frac{\partial y}{\partial \varphi} \\
\frac{\partial z}{\partial r} & \frac{\partial z}{\partial \theta} & \frac{\partial z}{\partial \varphi} \\
\end{vmatrix}
dr d\theta d\varphi 
```

```py3
X = r * sym.sin(t) * sym.cos(u)
Y = r * sym.sin(t) * sym.sin(u)
Z = r * sym.cos(t)
Jacobian = sym.Matrix([
    [sym.diff(X, r), sym.diff(X, t), sym.diff(X, u)],
    [sym.diff(Y, r), sym.diff(Y, t), sym.diff(Y, u)],
    [sym.diff(Z, r), sym.diff(Z, t), sym.diff(Z, u)]
])
2*sym.integrate(sym.integrate(sym.integrate(Jacobian.det().trigsimp(), (r, 0, 1)), (t, 0, Pi)), (u, 0, Pi/2))
```

$$\frac{2\pi}{3}$$

## 線積分
$$z=\frac{e^{-(x^2+y^2)}}{1+x^2+y^2}$$
を考える。

```py3
Z = sym.exp(-(x**2+y**2)) / (1+ x**2 + y**2)
plot3D(Z, (-2, 2, 0.1), (-2, 2, 0.1), 'surface')
```
![hiki5.png](https://qiita-image-store.s3.amazonaws.com/0/133885/3924517f-f41c-54bd-33e6-e45c031b60c3.png)
体積を求めると

```py3
X = r * sym.cos(t)
Y = r * sym.sin(t)
Z2 = Z.subs([(x, X), (y, Y)]).trigsimp()
Jacobian = sym.Matrix([[sym.diff(X, r), sym.diff(X, t)], [sym.diff(Y, r), sym.diff(Y, t)]])
result = sym.integrate(sym.integrate(Z2 * Jacobian.det().trigsimp(), (r, 0, oo)), (t, 0, 2*Pi))
sym.N(result, 30)
```
$1.87348049246219715751632002664$

半径1の円周上で線積分すると、

```py3
C = sym.Curve([sym.cos(t), sym.sin(t)], (t, 0, 2*Pi)) # 経路を媒介変数表示で与える
sym.line_integrate(Z, C, [x, y])
```
$$\frac{\pi}{e}$$

これを線積分の定義通りにやってみる。
媒介変数をそのまま代入して、
微分の自乗和の平方根を乗じて積分する。

```py3
X = sym.cos(t)
Y = sym.sin(t)
dX = sym.diff(X, t)
dY = sym.diff(Y, t)
sym.integrate(Z.subs([(x, X), (y, Y)])*sym.sqrt(dX**2 + dY**2), (t, 0, 2*Pi))
```
$$\frac{\pi}{e}$$