SymPy

SymPyの基本的な使い方

　SymPy とは Symbol (記号) による演算、すなわちコンピュータで代数演算を実行する機能を備えた計算機代数システム (CAS ; Computer Algebra System) です。

記号の定義

　x を記号として定義すると、x ** 2 + 1 のような式を x の多項式として扱うことができます。記号の定義には次の 2 通りの方法があります。

　複数の記号を定義することもできます。

　生成されたオブジェクトは sympy.core.symbol.Symbolクラスに属しますが、当サイトではこのクラスに属するオブジェクトを簡単に Symbol と記述します。

# sympyをインポート
import sympy

# Symbolクラスのインスタンスを生成（記号の定義）
sympy.var('x')

# 多項式yを定義
y = x ** 2 + 1

print(y)
print(type(x))

x**2 + 1
<class 'sympy.core.symbol.Symbol'>

記号に数値を代入

　記号に数値を代入するときは subs()メソッドを用います。

# sympyをインポート
import sympy

# Symbolクラスのインスタンスを生成
sympy.var('x')

# 多項式yを定義
y = x ** 2 + 1

# xに3を代入
y = y.subs(x, 3)

print(y)

10

　複数の記号に数値を代入するときは、(記号, 数値) のタプルを要素に持つリストを渡します。

# sympyをインポート
import sympy

# 記号 x, y を定義
sympy.var('x y')

# 多項式zを定義
z = x ** 2 + 2 * x * y + y ** 2

# xに3, yに5を代入
z = z.subs([(x, 3), (y, 5)])

print(z)

64

多項式の四則演算

　式同士で加減乗除を実行すると、多項式が返ります。

# sympyをインポート
import sympy

# 記号 a, b を定義
sympy.var('a b')

# 多項式f1を定義
f1 = a + 2 * b + 1
f2 = a + b + 5

# 多項式の加減乗除
g1 = f1 + f2
g2 = f1 - f2
g3 = f1 * f2
g4 = f1 / f2

print("f1 + f2 =", g1)
print("f1 - f2 =", g2)
print("f1 × f2 =", g3)
print("f1 / f2 =", g4)

f1 + f2 = 2*a + 3*b + 6
f1 - f2 = b - 4
f1 × f2 = (a + b + 5)*(a + 2*b + 1)
f1 / f2 = (a + 2*b + 1)/(a + b + 5)

分数型（有理数型）オブジェクト

　sympy.Rational() に分子と分母を渡すと、分数型（有理数型）インスタンスを生成することができます。正確なクラス名は sympy.core.numbers.Rational です。

# sympyをインポート
import sympy

# a = 3/5
a = sympy.Rational(3, 5)

# b = 1/7
b = sympy.Rational(1, 7)

# 記号を定義
sympy.var('c')

print("a =", a)
print("b =", b)
print("a + b + c =", a + b + c)

a = 3/5
b = 1/7a + b + c = c + 26/35

SymPy のバージョン確認とアップデート

　インストールされている SymPy のバージョンを確認したいときは __version__ 属性にアクセスしてください。

import sympy

# 現在の環境におけるSymPyのバージョンを確認する
print(sympy.__version__)

1.4

　SymPy を最新バージョンにアップデートする場合は、Anaconda プロンプトで以下のコマンドを打ちこみます。

> conda update sympy

　特定のバージョンをインストールする場合は以下のコマンドを実行してください。

> conda install sympy=1.2

　ただし、Anaconda 環境においては、単独パッケージのバージョン変更が他のライブラリとの依存関係を壊してしまう可能性があります。定期的に Anaconda 全体をアップデートすることをおすすめします。
　

検算にSymPyを使っています

　数学サイトの演習問題コーナーには BlogCat が作成したオリジナル問題が載っていたりしますが、万が一にも解答が間違っていたら大変なので（←何度も間違えているくせに）、ミスをなるべく減らすために、検算に Python を活用しています。一例として、演習問題 SQ-28 を題材に Python の実践コードを紹介します。SQ-28 では複素数の和を計算します。
　Python において、複素数は complex型として扱います。
　下のサンプルは $2k-1+i^k$ の $k=10$ ～ $12$ の和をとる計算です。

# sの初期値
s = 0

# 虚数iを定義
cpx = 1j

# a_k = 2k - 1 + i**k の k = 10, 11, 12 の和をとる
for k in range(10, 13):
    s += 2 * k - 1 + cpx ** k

print(s)

(63-1j)

　complex型は a + bj のような形式で定義します。
　上のコードでは虚数 $i$ を cpx = 1j と定義してあります。わざわざ変数を使う理由は、1j のような文字が式の中に入り込んでしまうと、

　　s += 2 * k - 1 + (1j) ** k

のようになって、目がちかちかしてして見えにくいからです。

　for文では range関数を使って 10, 11, 12 という数を順に生成しています。つい忘れがちですが、このような場合は range(10, 12) ではなく、range(10, 13) と記述します。

　次も複素数の和をとる計算です。先ほどより少し複雑で、
　
\[\sum_{M=1}^{12}\left\{M^2-\frac{(i-1)(i^M-1)}{2}\right\}\]
という計算をさせるコードです。

# sの初期値
s = 0

# 虚数iを定義
cpx = 1j

for k in range(1, 13):
    s += k ** 2 - (cpx - 1) * (cpx ** k - 1) / 2

print(s)

(644+6j)

　$k$ について $1$ から $12$ の和をとるので、for文のイテラブルオブジェクトは range(1, 13) です。このコードで虚数を変数に置き換えていないと、

　　s += k ** 2 - ((1j) - 1) * ((1j) ** k - 1) / 2

のようになって、目がちかちかちかっとなってしまいます。