【Python】sympyで数式の計算を自動化する

PythonのTopに戻る

 sympyとは

「sympy」とは、Pythonで記号計算を行うためのライブラリ。数値計算ではなく代数的に式の計算が可能であり、計算の検算などにも利用できる点が魅力的である。sympyは微分方程式や幾何、組み合わせ論などにも拡張されており、数学の諸分野を手広くカバーしている。詳細については公式サイトを確認して頂きたい。

例えば以下のような代数計算を自動で行ってくれる。

import sympy

x, y = sympy.symbols("x y")
f = (x+y)**3
print(sympy.expand(f))
# x**3 + 3*x**2*y + 3*x*y**2 + y**3

これより$$(x+y)^3=x^3 + 3 x^2 y + 3 x y^2 + y^3$$という関係式を瞬時に手に入れることができる。

なお、ネット上の情報を色々見てみると

# sympy.init_printing()

といった行が書かれているサンプルコードによく遭遇するが、これはJupyter Notebookで数式を表示する（レンダリングする）ための設定であり、一般のIDE上でsympyを利用する際に必要という訳ではない。

 シンボルの宣言

sympyにより文字を扱うには、まずシンボルの宣言を行う必要がある。先ほどのコードでも

x, y = sympy.symbols("x y")

という行が最初に書かれているが、これにより使用する文字をsympyモジュールに読み込ませている。このように複数の文字を一度に宣言するのが普通だが、添え字を含む文字を複数宣言する場合はこの方法では面倒である。そこで、愚直ではあるが例えば以下のような方法が思い付く。

import sympy

for i in range(0, 10):
    moji = 'x_{' + str(i) + '}'
    sympy.var(moji)
    moji_list.append(moji)
symbols = [sympy.symbols(v) for v in moji_list]
print(symbols)
# [x_{0}, x_{1}, x_{2}, x_{3}, x_{4}, x_{5}, x_{6}, x_{7}, x_{8}, x_{9}]

こうすることで、$x_{0}$から$x_{9}$までの10個の文字変数が一気に用意できる。しかしスライス記法を使えばfor文が不要となり、さらに簡単に書ける。

import sympy

a=sympy.symbols("a_{0:10}")
var = list(a)
print(var)
# [a_{0}, a_{1}, a_{2}, a_{3}, a_{4}, a_{5}, a_{6}, a_{7}, a_{8}, a_{9}]

波括弧は有っても無くても挙動上の違いは無いが、LaTeXコマンドとして利用する際に波括弧で添え字を括っておいた方が何かと便利である。

また、LaTex記法を使ってシンボルを宣言することもできる。このとき、raw文字列にしないとエスケープ文字と認識されるので注意する。

import sympy

theta, gamma = sympy.symbols(r"\theta \gamma")
f = sympy.sin(theta + gamma)
print(sympy.latex(f))

宣言された文字変数はsympyのデフォルトでは複素数として認識される。実数や整数、正の数、正整数の範囲で定義するには以下のようにする。

# 実数
a = sympy.symbols("a", real=True)
# 正の実数
b = sympy.symbols("b", positive=True)
# 整数
c = sympy.symbols("c", Integer=True)
# 正の整数
d = sympy.symbols("d", Integer=True, positive=True)
# 有理数
r = sympy.symbols("r", Rational=True)

数が属する集合を指定しておかないと求解の際に得られる結果が変わってくる。

 数式のLaTeXコマンド出力例

sympyの優れている点の一つが、数式のLaTeXコマンドを出力できる点である。例えば先ほどのコードにおいて出力するオブジェクトを sympy.latex(symbols) とすれば$$\left[ x_{0}, \ x_{1}, \ x_{2}, \ x_{3}, \ x_{4}, \ x_{5}, \ x_{6}, \ x_{7}, \ x_{8}, \ x_{9}\right]$$というLaTeX形式の表現が容易に得られる。

import sympy

for i in range(0, 10):
    moji = 'x_{' + str(i) + '}'
    sympy.var(moji)
    moji_list.append(moji)
symbols = [sympy.symbols(v) for v in moji_list]
print(sympy.latex(symbols))
# \left[ x_{0}, \ x_{1}, \ x_{2}, \ x_{3}, \ x_{4}, \ x_{5}, \ x_{6}, \ x_{7}, \ x_{8}, \ x_{9}\right]

計算結果のLaTeXコマンドを即座に入手したい場合は出力に sympy.latex() を作用させる。

import sympy

x, y = sympy.symbols("x y")
f = (x+y)**3
exp = sympy.expand(f)
fac = sympy.factor(exp)
print(sympy.latex(exp))
# x^{3} + 3 x^{2} y + 3 x y^{2} + y^{3}
print(sympy.latex(fac))
# \left(x + y\right)^{3}

極限計算は式に .limit(x, 0) などを作用させればよい。これは$\displaystyle \lim_{x \to 0}$と同じ意味。

import sympy

x = sympy.symbols("x")

f1 = sympy.sin(x)/x
print(sympy.latex(f1.limit(x, 0)))
# 1

f2 = (1-sympy.cos(x))/x**2
print(sympy.latex(f2.limit(x, 0)))
# \frac{1}{2}

定義式に対して .diff(x) を付ければ$x$で微分することも可能。この場合は偏微分になる。

import sympy

x, y = sympy.symbols("x y")
f = x * (x+y)**3
dif = f.diff(x)
exp = sympy.expand(f.diff(x))
fac = sympy.factor(exp)
print(sympy.latex(dif))
# 3 x \left(x + y\right)^{2} + \left(x + y\right)^{3}
print(sympy.latex(exp))
# 4 x^{3} + 9 x^{2} y + 6 x y^{2} + y^{3}
print(sympy.latex(fac))
# \left(x + y\right)^{2} \left(4 x + y\right)

これは$$f(x,y)=x(x+y)^3$$という関数を$x$について微分するプログラムである。これより、$$\begin{align} \dfrac{d}{dx}f(x,y) &=3 x \left(x + y\right)^{2} + \left(x + y\right)^{3} \\ &=\left(x + y\right)^{2} \left(4 x + y\right) \end{align}$$であり、展開すると $4 x^{3} + 9 x^{2} y + 6 x y^{2} + y^{3}$ となることが分かる。

また、sympyでは多くの関数について積分することも可能である。

import sympy

x = sympy.symbols("x")
for i in range(1, 6):
    f = sympy.sin(x)**i
    print("\["+sympy.latex(sympy.integrate(f, x)).replace("frac","dfrac")+"\]")

これは $\sin^n x$ を $n=1,2,…,5$ について積分した結果を出力する。\[- \cos{\left(x \right)}\] \[\dfrac{x}{2}-\dfrac{\sin{\left(x \right)} \cos{\left(x \right)}}{2}\] \[\dfrac{\cos^{3}{\left(x \right)}}{3}-\cos{\left(x \right)}\] \[\dfrac{3 x}{8}-\dfrac{\sin^{3}{\left(x \right)} \cos{\left(x \right)}}{4}-\dfrac{3 \sin{\left(x \right)} \cos{\left(x \right)}}{8}\] \[- \dfrac{\cos^{5}{\left(x \right)}}{5} + \dfrac{2 \cos^{3}{\left(x \right)}}{3}-\cos{\left(x \right)}\]

sympyでは線形代数も一通りカバーしている点も心強い。

import sympy

x, y = sympy.symbols("x y", real=True)
A = sympy.Matrix([[1,x], [y,1]])

print(sympy.latex(A**2))
# \left[\begin{matrix}x y + 1 & 2 x\\2 y & x y + 1\end{matrix}\right]

print(sympy.latex(sympy.simplify(A.inverse_CH())).replace("frac","dfrac"))
# \left[\begin{matrix}- \dfrac{1}{x y - 1} & \dfrac{x}{x y - 1}\\\dfrac{y}{x y - 1} & - \dfrac{1}{x y - 1}\end{matrix}\right]

これより$A=\left[\begin{matrix}1 & x \\ y & 1\end{matrix}\right]$ に対して、$$A^2=\left[\begin{matrix}x y + 1 & 2 x \\ 2 y & x y + 1\end{matrix}\right]$$と求められる。また、逆行列が$$\left[\begin{matrix}1 & x \\ y & 1\end{matrix}\right]^{-1}=\left[\begin{matrix}- \dfrac{1}{x y- 1} & \dfrac{x}{x y- 1}\\\dfrac{y}{x y- 1} &-\dfrac{1}{x y- 1}\end{matrix}\right]$$となることも分かる。行列計算はnumpyやscipyなどでも可能だが、sympyでは代数的に文字を処理できるので意外と有難い。

なお、逆行列の行に sympy.simplify を噛ませているが、これが無いと$$\scriptsize \left[\begin{matrix}\dfrac{- \dfrac{\left(x + y\right) \left(x-\dfrac{x + y}{y^{2} + 1}\right)}{\sqrt{y^{2} + 1} \left(x^{2}-\dfrac{\left(x + y\right)^{2}}{y^{2} + 1} + 1\right)} + \dfrac{1}{\sqrt{y^{2} + 1}}}{\sqrt{y^{2} + 1}} & \dfrac{\dfrac{y}{\sqrt{y^{2} + 1}}-\dfrac{\left(x + y\right) \left(- \dfrac{y \left(x + y\right)}{y^{2} + 1} + 1\right)}{\sqrt{y^{2} + 1} \left(x^{2}-\dfrac{\left(x + y\right)^{2}}{y^{2} + 1} + 1\right)}}{\sqrt{y^{2} + 1}}\\\dfrac{x-\dfrac{x + y}{y^{2} + 1}}{x^{2}-\dfrac{\left(x + y\right)^{2}}{y^{2} + 1} + 1} & \dfrac{- \dfrac{y \left(x + y\right)}{y^{2} + 1} + 1}{x^{2}-\dfrac{\left(x + y\right)^{2}}{y^{2} + 1} + 1}\end{matrix}\right]$$という約分される前の煩雑な式が出力される。

 その他の演算

その他の演算を紹介する。

・代入計算

import sympy

x, y, a = sympy.symbols("x y a", real=True)
A = sympy.Matrix([[x,y], [y,1]])

print(sympy.latex(A**2))

# y に 2a を代入
print(sympy.latex(A.subs([(y, 2*a)])**2))

$$A^2=\left[\begin{matrix}x^{2} + y^{2} & x y + y \\ x y + y & y^{2} + 1\end{matrix}\right]$$に対して $y=2a$ を代入すると$$A^2=\left[\begin{matrix}4 a^{2} + x^{2} & 2 a x + 2 a \\ 2 a x + 2 a & 4 a^{2} + 1\end{matrix}\right]$$を得る。

・方程式を解く

import sympy

x, y = sympy.symbols("x y", real=True)
eq = sympy.Eq(2*x - 3*y, 5)
print(sympy.solve(eq, x))
# [3*y/2 + 5/2]

これは $2x-3y=5$ という方程式の解を求めている。

import sympy

x, y = sympy.symbols("x y", real=True)
eq = sympy.Eq(2*x**2 - 3*y, 5)
print(sympy.solve(eq, x))
# [-sqrt(6*y + 10)/2, sqrt(6*y + 10)/2]

これは $2x^2-3y=5$ という方程式の解を求めている。一般に、特殊な形の方程式を除いて、5次以上の方程式の解はsympyでは得られない。

連立方程式は次のように解くことができる。

import sympy

x, y, z = sympy.symbols("x y z", Integer=True)
print(sympy.linsolve([x + y - 2*z - 1, x + y + 2*z - 3 ], (x, y, z)))
# FiniteSet((2 - y, y, 1/2))
print(sympy.latex(sympy.linsolve([x + y - 2*z - 1, x + y + 2*z - 3 ], (x, y, z))))
# \left\{\left( 2 - y, \  y, \  \frac{1}{2}\right)\right\}

これは$$\begin{cases} x+y-2z=1 \\ x+y+2z=3 \end{cases}$$という3元1次連立方程式の求解に相当し、解の集合として$\left\{\left( 2-y, \ y, \ \dfrac{1}{2}\right)\right\}$が得られる。

・総和＆総乗計算

import sympy

i, N = sympy.symbols("i N", Integer=True)
print(sympy.summation(i**2, (i, 1, N)))
# N**3/3 + N**2/2 + N/6
print(sympy.latex(sympy.summation(i**2, (i, 1, N))).replace("frac","dfrac"))
# \dfrac{N^{3}}{3} + \dfrac{N^{2}}{2} + \dfrac{N}{6}

これは$$\displaystyle \sum_{i=1}^{N} i^2=\dfrac{N^{3}}{3} + \dfrac{N^{2}}{2} + \dfrac{N}{6}$$を求めるコードである。$\displaystyle \sum_{i=1}^{N} \sqrt{i}$ のように閉じた式で表せない和の場合は Sum(sqrt(i), (i, 1, N)) という値を返す。

for文と組み合わせれば総和公式のリストが得られる。

$\small \dfrac{1}{2}n \left(n + 1\right)$
$\small \dfrac{1}{6}n \left(n + 1\right) \left(2 n + 1\right)$
$\small \dfrac{1}{4}n^{2} \left(n + 1\right)^{2}$
$\small \dfrac{1}{30}n \left(n + 1\right) \left(2 n + 1\right) \left(3 n^{2} + 3 n-1\right)$
$\small \dfrac{1}{12}n^{2} \left(n + 1\right)^{2} \left(2 n^{2} + 2 n-1\right)$
$\small \dfrac{1}{42}n \left(n + 1\right) \left(2 n + 1\right) \left(3 n^{4} + 6 n^{3}-3 n + 1\right)$
$\small \dfrac{1}{24}n^{2} \left(n + 1\right)^{2} \left(3 n^{4} + 6 n^{3}-n^{2}-4 n + 2\right)$
$\small \dfrac{1}{90}n \left(n + 1\right) \left(2 n + 1\right) \left(5 n^{6} + 15 n^{5} + 5 n^{4}-15 n^{3}-n^{2} + 9 n-3\right)$
$\small \dfrac{1}{20}n^{2} \left(n + 1\right)^{2} \left(n^{2} + n-1\right) \left(2 n^{4} + 4 n^{3}-n^{2}-3 n + 3\right)$
$\small \dfrac{1}{66}n \left(n + 1\right) \left(2 n + 1\right) \left(n^{2} + n-1\right) \left(3 n^{6} + 9 n^{5} + 2 n^{4}-11 n^{3} + 3 n^{2} + 10 n-5\right)$

また、等比級数でも余程煩雑な式でない限り総和計算は上手くいく。

import sympy

k, N = sympy.symbols("k N", Integer=True)
print(sympy.factor(sympy.summation(3**k, (k, 1, N) )))
# 3*(3**N - 1)/2

$$\displaystyle \sum_{k=1}^{N} 3^k = \dfrac{3 \left(3^{N}-1\right)}{2}$$

print(sympy.summation(sympy.sin(sympy.Rational(2,3) *n*sympy.pi)/(10**n), (n, 1, sympy.oo)))

なお、総乗計算には sympy.product 関数を用いる。

import sympy

k, N = sympy.symbols("k N", Integer=True)
print(sympy.factor(sympy.product(3**k, (k, 1, N) )))
# 3**(N**2/2 + N/2)

$$\displaystyle \prod_{k=1}^{N} 3^k = 3^{N(N + 1)/2}$$と正しく求められている。

・近似値の表示

import sympy

print(sympy.pi)
# pi
print(float(sympy.pi))
# 3.141592653589793

print(float(sympy.exp(1)))
print(float(sympy.E))
# 2.718281828459045

float関数を使う。自然対数の底（ネイピア数）$e$は sympy.E で定義され、虚数単位は sympy.I で定義される。以下のようにしてオイラーの等式 $e^{i \pi}=-1$ の成立が確かめられる。

import sympy

print(sympy.E**(sympy.I * sympy.pi))
# -1

PythonのTopに戻る