数式操作で便利なSymPyの入門

数式表現や操作・可視化などで便利なSymPyライブラリについて入門していくための記事です。

※記事執筆者は理系出身ではなく数学や用語などで色々粗い点等はご容赦ください。

# SymPyでできること

以下一例です。後々の節で順番に触れていきます。

- Python上での数式表示が楽になります。
- Python上での数式表現や演算がやりやすくなります。
- ダイレクトに数式を使った計算を行うことができます。
- 数式の解を求めることができます。
- 数式のプロットを行うことができます（内部ではmatplotlibが使われます）。

# 使うもの

- Windows 10
- Python 3.9.0
- sympy==1.9
- Jupyter（VS Code上のものを利用しているため背景などが黒くなっています）

# SymPyのインストール

pipで対応することができます。

```
$ pip install sympy
```

※今回は記事執筆時点で最新の1.9のバージョンを使っているため合わせる場合には以下のようなコマンドをお使いください。

```
$ pip install sympy==1.9
```

# シンボルの定義

$x$や$y$、$a$といったようなアルファベットのシンボルを定義するにはSymPyのシンボルクラスを使います。第一引数にはシンボルに設定するアルファベットを指定します。基本的に変数名と同じ名前を設定します。

例えば$x$というシンボルを定義したい場合には以下のように書きます。

```py
from sympy import Symbol

x: Symbol = Symbol('x')
```

生成したシンボルをそのままJupyter上で表示すると数式表示のアウトプットになります。少し分かりづらいですがアウトプットでセンタリングされる形で$x$と表示されています。

```py
x
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/228778/fd32048d-8fec-870b-4c24-d85203d6b2d8.png)

もしくはsymbols関数でまとめてシンボルを定義することもできます。引数にコンマ区切りで各シンボル名を指定します。ただしこちらの場合は現状Pylance（Pyright）などでの型チェックに引っかかります。動作上は問題無く動作します。

```py
from sympy import symbols
x, y = symbols('x, y')
```

# シンボルに対する演算

基本的な演算処理を行う各オペレーターは直接シンボルに対して行うことができます。例えば以下のように2乗や乗算、加算などを行うことができます。

```py
expr = x ** 2 + 2 * x + 1
```

演算結果はJupyter上で表示してみると数式としてちゃんと反映されていることを確認できます。

```py
expr
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/228778/d175cd0d-72d7-f9d0-e13a-110621f793e0.png)

括弧などを付けてもちゃんと数式に反映されます。

```py
expr = (x + 1) * (x + 2)
expr
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/228778/9b012d56-352e-3deb-a198-38f5d938bc24.png)

# シンボル名の取得や表示

シンボル名を文字列として取得したり表示したい場合にはいくつか方法があります。

1つ目としてはname属性にアクセスする方法です。

```py
name: str = x.name
print('name:', name, '\ntype:', type(name))
```

```
name: x
type: <class 'str'>
```

2つ目としてstr型にキャストする方法も使えます。こちらもそのままシンボル名の文字列になります。

```py
name: str = str(x)
print('name:', name, '\ntype:', type(name))
```

```
name: x
type: <class 'str'>
```

また、print関数にそのまま通してもシンボル名が表示されます。

```py
print(x)
```

```
x
```

# 数式の展開

expand関数で数式の展開を行うことができます。

例として$(x + y)^2$という式を展開してみます。まずは対象の式を用意します。

```py
x = Symbol('x')
y = Symbol('y')
expr = (x + y) ** 2
expr
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/228778/c6fc1159-8727-1157-00e0-74c4090e659d.png)

この式をexpand関数を使って展開してみます。

```py
from sympy import expand

expr = expand(expr)
expr
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/228778/0da34c3e-8154-c3bf-cbcb-5eee5a0e98ab.png)

$x^2 + 2xy + y^2$というように展開してくれました。

# 数式の因数分解

今度は展開とは逆に因数分解してみます。factor関数で扱うことができます。

先ほどの$x^2 + 2xy + y^2$という式をまずは用意します。

```py
expr = x ** 2 + 2 * x * y + y ** 2
expr
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/228778/f90ad7ce-b64b-cb45-f499-791dcba54ce3.png)

factor関数を通してみます。

```py
from sympy import factor

expr = factor(expr)
expr
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/228778/94d77ad6-a12d-0769-0635-0336991e8c3d.png)

展開前の式に戻すことができました。

# 式の自動の計算の挙動

SymPyではシンプルな表現部分は自動で計算され統合などがされます。例えば以下のように$x + x$というシンプルな式であれば$2x$という表示になります。

```py
x + x
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/228778/c8ee3045-ead8-30b8-4c34-6195533a4215.png)

一方で式が複雑になってくると計算を行わずにそのまま扱った方が良いケースなども多くなってきます。そういった計算はそのまま残る形となります（全て展開されるわけではありません）。

例えば以下の式では$x + x$の部分が$2x$、$2 * x$の部分が$2x$という形になっていますが、一方で結果は$2x(2x + 1)$となり全て計算されているわけではありません。

```py
(x + x) * (2 * x + 1)
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/228778/d79a3b58-d119-abdc-3c68-90017ea5ddd6.png)


# 式の要素の表示順

デフォルトではSymPyでの表示は指数の値が大きい方が左側に表示されます。例えば以下のようにPython上で書くとPython上の記述と反対の順番で式が表示されます。

```py
5 + x + x ** 2 + x ** 3
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/228778/1cf82cc6-c0e2-5481-e9e4-f7acf860a002.png)

通常は式の表示としてはこれが自然ですが、何らかの理由で逆順に表示したい場合はinit_printing関数で`order='rev-lex'`と指定します。rev-lexはreverse lexicographical order （辞書的に逆順）の略となります。この辺は詳しくは触れませんが必要な方は以下のWikipediaの記事などをご確認ください。

https://en.wikipedia.org/wiki/Lexicographic_order

```py
from sympy import init_printing
init_printing(order='rev-lex')
5 + x + x ** 2 + x ** 3
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/228778/2d0f96e8-5264-3417-ea04-a79fc1d429d0.png)

# ターミナルやコマンドプロンプト上などでも見栄え良く表示したい場合には・・・

Jupyter上などで扱っているとTeXで表示してくれるため特に意識しなくとも良い感じに表示してくれます。一方でターミナルやコマンドプロンプト上での表示が必要な場合にはこの表示は使えません。このようなケースでもある程度読みやすくするためのSymPyのpprint関数が用意されています。

例えば以下のような数式で考えてみます。

```py
expr = 1 + x ** 2 / 2 + x ** 3 / 3 + x ** 4 / 4 + x ** 5 / 5
```

Jupyter上で表示してみると以下のように表示されます。

```py
expr
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/228778/1c35f5e5-d5d5-04d1-b618-e1432a2e312c.png)

ターミナル上で表示したり、もしくはJupyter上などでprintで表示すると以下のようにテキストで表示されます。ちょっと読みやすい・・・とは言い難い表示となります。

```py
print(expr)
```

```
x**5/5 + x**4/4 + x**3/3 + x**2/2 + 1
```

こういった場合にはSymPyのpprint（pretty-print）を使うと改行やスペースなどが良い感じに入って先ほどのものよりも読みやすい表示になります（指数部分などが小さくなったりは流石に無理ですが・・・）。

```py
from sympy import pprint

pprint(expr)
```

```
     2    3    4    5
    x    x    x    x 
1 + ── + ── + ── + ──
    2    3    4    5 
```

# 各シンボルに代入を行う

※init_printingの設定をリセットするため一旦Jupyterのカーネルを再起動しています。

式表現のオブジェクトはsubsメソッドを持っているため、そちらに値を指定することで各シンボルに代入を行うことができます（substitutionで代入という意味になります）。

例えば以下のような式で考えてみます。

```py
from sympy import Symbol

x = Symbol('x')
expr = 3 * x ** 2 + 5 * x + 2
expr
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/228778/028b3d55-b739-6768-a08e-88a06fa81572.png)

この式に対してsubsメソッドでx = 6という値を代入してみましょう。辞書の形でキーに対象のシンボル、値に代入したい値を指定して使います。

```py
expr.subs({x: 6})
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/228778/3f83f847-b294-a801-e4cf-fc242cfea9b4.png)

3 * 6 ** 2 + 5 * 6 + 2で140という値を得ることができました。

この代入処理は別のシンボルを含めることもできます。例えばx = y + 5という条件だった場合、そのまま辞書の値に`y + 5`という値を指定することができます。

```py
y = Symbol('y')
expr.subs({x: y + 5})
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/228778/5f4028d3-9b12-dff7-f6be-c2a1b1791235.png)

# 文字列からSymPyのオブジェクトへ変換する

DBへ保存したりユーザーからの入力値を受け付ける場合などには値がSymPyの数式のオブジェクトではなく数式表現用の文字列になります。これらの文字列をSymPyのオブジェクトへ変換したり戻したりしたい場合にはsympifyという関数が用意されておりそちらを使うことで実現できます。

```py
from sympy import sympify

expr = sympify('3 * x ** 2 + 5 * x + 2')
expr
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/228778/5a99c260-66f4-1e18-9bad-8c759e506e12.png)

文字列から変換できない場合（不正な数式表現の文字列の場合）などにはエラーになります。

例えば以下のように文字列上では`3 * x`と`*`を挟まないといけない箇所で`3x`などとなっていた場合などにはエラーになります。

```py
expr = sympify('3x')
```

```
SympifyError: Sympify of expression 'could not parse '3x'' failed, because of exception being raised:
SyntaxError: invalid syntax (<string>, line 1)
```

# 式の解を求める

式でどんな値を設定すると0になるか・・・といった解を求めることがsolve関数を使って行うことができます。

例えば$x^2 - 4x + 3 = 0$という二次方程式の解を考えてみます。

まずは数式を用意します。

```py
expr = x ** 2 - 4 * x + 3
expr
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/228778/49bec718-eded-02a7-b8c4-c19057eb276e.png)

そちらを使ってsolve関数を通してみます。

```py
from sympy import solve
solve(expr)
```

```
[1, 3]
```

1と3という解が得られました。なお今回使った2次方程式では解が2つですが、解が1つの式の場合でも結果はリストで返ってきます。

試しに得られた1と3の解をsubsメソッドで代入して結果が0になることを確認してみます。

```py
expr.subs({x: 1})
```

```py
expr.subs({x: 3})
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/228778/901c25bc-263c-a906-2e55-ed21510c81cd.png)

それぞれの解を代入すると結果が0になることを確認することができました。

# 数式をプロットする

SymPyでは内部でmatplotlibを使う形でプロットするインターフェイスが用意されているため、数式をダイレクトにプロットすることができます。

例として$y = 3x + 5$という一次関数の式を使ってみます。

```py
expr = 3 * x + 5
expr
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/228778/b6c0948d-5bc2-1887-fa42-832bec0007ea.png)

プロットのインターフェイスは`sympy.plotting`パッケージに`plot`関数が含まれているのでそちらを使います。

```py
from sympy.plotting import plot
plot(expr)
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/228778/1bd2cf99-0a32-ff22-9a66-44ee5850a1a8.png)

可視化までがお手軽ですね。

## 横軸の範囲を調整する

plot関数は前節で可視化された結果の通り、デフォルトでは横軸（前節ではxの値）は-10～10の範囲で設定されます。

この範囲を変えたい場合にはplot関数の引数に2件の引数を指定します。第一引数には対象の数式、第二引数には対象のシンボル・横軸最小値・横軸最大値の3つの値を格納したタプルを指定します。

たとえば-5～5の範囲にしたければ以下のように書きます。

```py
from sympy.plotting import plot
plot(expr, (x, -5, 5))
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/228778/874c320c-455a-80c1-eb61-3089037d06cf.png)

## タイトルや軸のラベルなどを設定する

plot関数にはtitle, xlabel, ylabelといったmatplotlibでお馴染みの引数があるため、それらをキーワード引数で指定することもできます。

```py
from sympy.plotting import plot
plot(expr, title='Linear function of 3x + 5', xlabel='x value')
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/228778/e3081927-9860-f00f-cf0e-ff7d96ce02cf.png)

## 複数の数式を同時に表示する

複数の数式を同時にプロットしたい場合には第一引数に1つ目の数式のオブジェクト、第二引数に2つ目の数式のオブジェクト・・・と順番に指定していくことで対応できます。

```py
from sympy.plotting import plot
expr1 = 3 * x + 5
expr2 = 2 * x + 10
plot(expr1, expr2)
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/228778/248c19fe-c42a-35ee-4d8d-dcc19555be42.png)

# 参考文献・参考サイトまとめ

https://www.amazon.com/Doing-Math-Python-Programming-Statistics/dp/1593276400/

https://en.wikipedia.org/wiki/Lexicographic_order

https://yorikuwa.com/m1300/

https://ist.ksc.kwansei.ac.jp/~nishitani/Maple/JtoE.html