Pythonで数学の勉強：代数基礎(数式と関数)

# 数式と関数
## 目的など

* よくわからなかった大学の数学をPythonを使って再勉強する。(今ならMathematica並のことができるぜ)
* 数学自体が目的ではなく現実の問題を解くために数学の知識を利用したい。
* なるべく例題を入れつつ、解答は例題に特化したものではなく汎用的なものにする。
* どこがsympyの領域なのかを明確にすべくfrom sympy import *とはしない。

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

# 使用する変数の定義(小文字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')
```

## 式の定義

```py3
F1 = x * x + 2 * x + 1 # 普通にシンボルを演算するだけ
F2 = x ** 5 + x + 1 # 高次式
F3 = (x + 2) * (x - 1) / ((x - 2) * (x - 2)) # 有理関数
F4 = sym.cos(Pi * x) # 三角関数
F5 = sym.exp(x) / sym.Pow(x, 10) # 指数関数とべき乗
F6 = sym.log(x * x + 1) # 対数関数
F7 = (x + y) ** 3 # 2変数関数
F8 = (x + y + z / 2) ** 4 # 多変数関数
```

$x^{2} + 2 x + 1$
$x^{5} + x + 1$
$\\frac{\\left(x - 1\\right) \\left(x + 2\\right)}{\\left(x - 2\\right)^{2}}$
$\\cos{\\left (\\pi x \\right )}$
$\\frac{e^{x}}{x^{10}}$
$\\log{\\left (x^{2} + 1 \\right )}$
$\\left(x + y\\right)^{3}$
$\\left(x + y + \\frac{z}{2}\\right)^{4}$

## 式の展開

```py3
sym.expand(F8)
```
$x^{4} + 4 x^{3} y + 2 x^{3} z + 6 x^{2} y^{2} + 6 x^{2} y z + \\frac{3 x^{2}}{2} z^{2} + 4 x y^{3} + 6 x y^{2} z + 3 x y z^{2} + \\frac{x z^{3}}{2} + y^{4} + 2 y^{3} z + \\frac{3 y^{2}}{2} z^{2} + \\frac{y z^{3}}{2} + \\frac{z^{4}}{16}$

## 因数分解

```py3
sym.factor(F1)
```
$(x+1)^2$


## 式への代入

```py3
F4.subs([(x, sym.Rational(1, 4))])
F6.subs([(x, 5)])
F8.subs([(x, 1), (y, -1), (z, 3)])
```

$\\frac{\\sqrt{2}}{2}$
$\log (26)$
$\\frac{81}{16}$

## 式の簡略化
* powsimp (指数の簡略化), 
* trigsimp (三角関数を含む式の簡略化), 
* logcombine (対数関数を含む式の簡略化)
* radsimp (複素数を含む式の簡略化)

```py3
sym.radsimp((2+3*I)/(3-4*I))
```
$\frac{1}{25}(-6+17i)$

## 浮動小数への変換

```py3
print(sym.N(sym.sin(1), 30)) # 桁数を指定
print(sym.N(Pi, 50))
print(float(sym.N(Pi, 50))) # Python標準の形式に戻す
```
```text
0.841470984807896506652502321630
3.1415926535897932384626433832795028841971693993751
3.141592653589793
```

## LaTex式の出力

```py3
(F1*F2*F3+F4)/(F5*F6/F7+F8)
```
![siki3.png](https://qiita-image-store.s3.amazonaws.com/0/133885/e1d85487-76c4-2329-5178-871673b55f94.png)

```py3
sym.latex((F1*F2*F3+F4)/(F5*F6/F7+F8))
```
```text
\\frac{\\cos{\\left (\\pi x \\right )} + \\frac{1}{\\left(x - 2\\right)^{2}} \\left(x - 1\\right) \\left(x + 2\\right) \\left(x^{2} + 2 x + 1\\right) \\left(x^{5} + x + 1\\right)}{\\left(x + y + \\frac{z}{2}\\right)^{4} + \\frac{e^{x} \\log{\\left (x^{2} + 1 \\right )}}{x^{10} \\left(x + y\\right)^{3}}}
```

論文書くのに便利機能すぎる。

## 特定の変数について解く

```py3
F9 = x * x + sym.sin(y + z ** 3)
sym.solve(F9, z)
```
![siki4.png](https://qiita-image-store.s3.amazonaws.com/0/133885/011fbc52-4f6d-717e-4f30-8b56f3d5261a.png)

$F_9=0$を$z$について解いている。
まず$\sin$の逆関数で$\pm$、さらに三乗根で3倍されて計6個の解がある。

## 方程式を解く

```py3
sym.solve(F2, x) # F2 = x ** 5 + x + 1
```
![siki5.png](https://qiita-image-store.s3.amazonaws.com/0/133885/0a6aa980-c35c-62c9-0e94-b2cd54c3535a.png)

解けるんかい。
解を見ると1の3乗根2つがうまく打ち消しあって
これで割れば残りは3次式というわけか。
ちなみに$x^5 + 2x + 1$だと解析的な解は求まりません。

## 連立方程式を解く

```py3
F10 = 2*x + 2*y + 6*z + 3
F11 =   x +   y - 2*z - 7
F12 = 5*x - 9*y + 2*z - 15
sym.solve([F10, F11, F12], [x, y, z])
```
![siki6.png](https://qiita-image-store.s3.amazonaws.com/0/133885/5d354733-16dc-34d0-8908-a990d51da211.png)

リストを渡すだけ。簡単。

## 演習
1. $x=\frac{2+5i}{3-i}, y=\frac{i-2}{3i+1}$ のとき次の式の値を計算せよ。
$$\frac{1}{x^2 + 3xy - y^2}$$
2. 次の式を展開せよ。
$$(a+b+c)^2-(b+c-a)^2+(c+a-b)^2-(a+b-c)^2$$
3. 次の式を因数分解せよ。
$$a^2(y-b)-(a-b)xy+ab(x+y)+b^2(x-a)$$
4. 次の式をPythonの通常の方法とSympyを使った方法とで比較せよ。
$$\sum_{k=0}^{10000} \sin(\frac{k}{10000})$$

5. 次の方程式の解を求めよ。
$$\sin^2{x}-\cos^2{x}=0$$

6. 次の連立方程式を解け。
${x^2+y^2=1},\\  
{y=x+\frac{1}{2}}$

```py3
X = (2+5*I)/(3-I)
Y = (I-2)/(3*I+1)
F = 1/(x*x + 3*x*y - y*y)
sym.expand(sym.radsimp(F.subs([(x, X), (y, Y)])))
```
$$- \frac{594}{3613} - \frac{92 i}{3613}$$

```py3
sym.expand((a+b+c)**2-(b+c-a)**2+(c+a-b)**2-(a+b-c)**2)
```

$$8ac$$

```py3
sym.factor(a**2 *(y-b)-(a+b)*x*y+a*b*(x+y)+b**2 *(x-a))
```

$$−(−a+x)(a+b)(−b+y)$$

```py3
import math
print(math.fsum([math.sin(K/10000.0) for K in range(10000)])) # Python標準
```
4596.556201995385

```py3
# 10秒ぐらいかかります
r = 0
for K in range(10000):
    r += sym.N(sym.sin(sym.Rational(K, 10000)), 50) # なるべく内側で浮動小数にする
r
```
4596.5562019953847593333362996376237160871901535447

Pythonのfsumが思った以上に優秀だった。

```py3
F = (sym.sin(x) * sym.sin(x) - sym.cos(x) * sym.cos(x))
sym.solve(F, x)
```
$$\\left [ - \\frac{3 \\pi}{4}, \\quad - \\frac{\\pi}{4}, \\quad \\frac{\\pi}{4}, \\quad \\frac{3 \\pi}{4}\\right ]$$

```py3
X = x * x + y * y - 1
Y = x - y + sym.Rational(1, 2)
sym.solve([X, Y], [x, y])
```

$$\\left [ \\left ( - \\frac{1}{4} + \\frac{\\sqrt{7}}{4}, \\quad \\frac{1}{4} + \\frac{\\sqrt{7}}{4}\\right ), \\quad \\left ( - \\frac{\\sqrt{7}}{4} - \\frac{1}{4}, \\quad - \\frac{\\sqrt{7}}{4} + \\frac{1}{4}\\right )\\right ]$$

 円と直線の交点ですね。