SymPy完全入門：Pythonで記号計算・数式処理を自在に操る

はじめに：SymPyとは何か？

SymPyは、Pythonプログラミング言語で書かれた、記号計算（シンボリック計算）のためのオープンソースライブラリです。数値計算ライブラリであるNumPyやSciPyが具体的な数値を用いて計算を行うのに対し、SymPyは数学の式を文字や記号のまま扱います。これにより、代数的な操作、方程式の解析解の導出、微積分、行列演算などを記号的に行うことができます。

MathematicaやMapleのような高機能な商用数式処理システム（CAS: Computer Algebra System）の代替となることを目指しつつも、コードのシンプルさと拡張性を重視して開発されています。完全にPythonで書かれており、基本的な機能は外部ライブラリ（mpmathを除く）に依存しないため、導入が容易であるという特徴も持っています。

科学技術計算、教育、研究開発など、幅広い分野でその強力な機能が活用されています。

基本的な使い方

インストール

SymPyのインストールは、Pythonのパッケージマネージャであるpipを使うのが最も簡単です。ターミナルやコマンドプロンプトで以下のコマンドを実行します。

pip install sympy

Anacondaを使用している場合は、condaコマンドでもインストールできます。

conda install sympy

シンボル（記号）の定義

SymPyで数式を扱うには、まず変数として使用する記号（シンボル）を定義する必要があります。これはsympy.Symbolまたはsympy.symbols関数を使用します。

import sympy

# シンボル 'x' を定義
x = sympy.Symbol('x')

# 複数のシンボル 'y', 'z' を同時に定義
y, z = sympy.symbols('y z')

# シンボルに関する情報を付加することも可能 (例: 実数, 正の数)
a = sympy.Symbol('a', real=True)
b = sympy.Symbol('b', positive=True)

print(x)
print(y, z)
print(a.is_real, b.is_positive)

定義したシンボルを使って、Pythonの演算子で数式を組み立てることができます。

expr = x**2 + 2*y + z
print(expr)  # 出力: x**2 + 2*y + z

注意: Pythonの変数名（例: x）とSymPyのシンボル名（例: 'x'）は別物ですが、通常は同じ名前を対応付けて使うのが分かりやすいです。

基本的な演算と数式表現

定義したシンボルとPythonの標準的な演算子（+, -, *, /, **）を使って、数式を表現できます。SymPyはこれらの式を記号的に保持します。

expr1 = (x + y)**2
expr2 = sympy.sin(x) + sympy.cos(y)

print(expr1)  # 出力: (x + y)**2
print(expr2)  # 出力: sin(x) + cos(y)

円周率 π や自然対数の底 e、無限大 ∞ などの数学定数も用意されています。

print(sympy.pi)      # 出力: pi
print(sympy.E)       # 出力: E (自然対数の底)
print(sympy.oo)      # 出力: oo (無限大)
print(sympy.I)       # 出力: I (虚数単位)

# 定数を含む計算
area = sympy.pi * (x/2)**2
print(area)          # 出力: pi*x**2/4

式の評価と代入

数式中のシンボルに具体的な値を代入するには、subs()メソッドを使用します。

expr = x**2 + 3*x + 2

# x に 1 を代入
result1 = expr.subs(x, 1)
print(result1)  # 出力: 6

# x に y+1 を代入 (シンボルも代入可能)
result2 = expr.subs(x, y + 1)
print(result2)  # 出力: (y + 1)**2 + 3*(y + 1) + 2

# 複数の値を一度に代入 (辞書またはタプルのリストを使用)
expr_multi = x**2 + y**2
result3 = expr_multi.subs({x: 1, y: 2})
print(result3)  # 出力: 5

result4 = expr_multi.subs([(x, 1), (y, 2)])
print(result4)  # 出力: 5

シンボリックな表現を数値（浮動小数点数）として評価したい場合は、evalf()メソッドを使用します。精度を指定することも可能です。

pi_val = sympy.pi.evalf()
print(pi_val)  # 出力: 3.14159265358979 (デフォルトの精度)

pi_val_high = sympy.pi.evalf(50) # 50桁の精度で評価
print(pi_val_high) # 出力: 3.1415926535897932384626433832795028841971693993751

N()関数もevalf()と同様の機能を提供します。

sqrt2_approx = sympy.N(sympy.sqrt(2), 20) # 20桁の精度で評価
print(sqrt2_approx) # 出力: 1.4142135623730950488

出力形式の調整 (Pretty Printing)

Jupyter NotebookやIPython環境では、sympy.init_printing()を実行すると、数式がLaTeXのようなきれいな形式で表示されるようになります。

sympy.init_printing(use_unicode=True) # Unicode文字を使って表示を綺麗にする

expr = sympy.Integral(sympy.sqrt(1/x), x)
expr # Jupyter Notebook/IPython上できれいに表示される

特定の式をLaTeXコードとして取得したい場合は、sympy.latex()関数を使用します。

latex_code = sympy.latex(expr)
print(latex_code) # 出力: \int \sqrt{\frac{1}{x}}\, dx

主要な機能

SymPyは非常に多くの機能を提供していますが、ここでは特に重要なものをいくつか紹介します。

代数演算

• 展開 (Expand): sympy.expand() 関数で多項式や式を展開します。
• 因数分解 (Factor): sympy.factor() 関数で式を因数分解します。
• 簡略化 (Simplify): sympy.simplify() 関数で式をできるだけ簡単な形にします。これは様々なヒューリスティックを用いて式の簡略化を試みます。特定の種類の簡略化（三角関数、指数・対数など）に特化した関数もあります（例: trigsimp, powsimp, logcombine）。
• 整理 (Collect): sympy.collect() 関数で特定の変数について式を整理します。
• 部分分数分解 (Apart): sympy.apart() 関数で有理式を部分分数に分解します。

sympy.init_printing(use_unicode=True)
x, y = sympy.symbols('x y')

# 展開
expr1 = (x + y)**3
expanded_expr1 = sympy.expand(expr1)
print("展開:")
sympy.pprint(expanded_expr1) # pprint は整形して表示する関数
# 出力:
# 展開:
# 3      2        2      3
# x  + 3⋅x ⋅y + 3⋅x⋅y  + y

# 因数分解
expr2 = x**3 + 3*x**2*y + 3*x*y**2 + y**3
factored_expr2 = sympy.factor(expr2)
print("\n因数分解:")
sympy.pprint(factored_expr2)
# 出力:
# 因数分解:
#        3
# (x + y)

# 簡略化
expr3 = sympy.sin(x)**2 + sympy.cos(x)**2
simplified_expr3 = sympy.simplify(expr3)
print("\n簡略化:")
sympy.pprint(simplified_expr3)
# 出力:
# 簡略化:
# 1

# 特定の変数について整理
expr4 = x*y + x - 3 + 2*x**2 - z*x**2 + x**3
collected_expr4 = sympy.collect(expr4, x)
print("\nxについて整理:")
sympy.pprint(collected_expr4)
# 出力:
# xについて整理:
#  3      2              2
# x  + (-z + 2)⋅x  + x⋅y + x - 3

# 部分分数分解
expr5 = 1 / ((x + 1) * (x + 2))
apart_expr5 = sympy.apart(expr5, x)
print("\n部分分数分解:")
sympy.pprint(apart_expr5)
# 出力:
# 部分分数分解:
#    1       1
# ───── - ─────
# x + 1   x + 2

微積分

• 極限 (Limit): sympy.limit(式, 変数, 近づける値, 方向) で極限値を計算します。
• 微分 (Derivative): sympy.diff(式, 変数, 回数) で微分を計算します。多変数関数の偏微分も可能です。
• 積分 (Integral): sympy.integrate(式, 変数) で不定積分を、sympy.integrate(式, (変数, 下端, 上端)) で定積分を計算します。多重積分も可能です。
• テイラー級数展開 (Series Expansion): 式.series(変数, 点, 次数) で指定した点の周りでのテイラー級数（ローラン級数）展開を求めます。

sympy.init_printing(use_unicode=True)
x, y, h = sympy.symbols('x y h')

# 極限: lim (sin(x)/x) as x -> 0
limit_val = sympy.limit(sympy.sin(x)/x, x, 0)
print("極限:")
sympy.pprint(limit_val)
# 出力:
# 極限:
# 1

# 微分: d/dx (x^3 * sin(x))
diff_expr = sympy.diff(x**3 * sympy.sin(x), x)
print("\n微分:")
sympy.pprint(diff_expr)
# 出力:
# 微分:
#  3           2
# x ⋅cos(x) + 3⋅x ⋅sin(x)

# 偏微分: ∂^2/∂x∂y (exp(x*y))
partial_diff_expr = sympy.diff(sympy.exp(x*y), x, 1, y, 1) # xで1回、yで1回偏微分
print("\n偏微分:")
sympy.pprint(partial_diff_expr)
# 出力:
# 偏微分:
#     x⋅y      x⋅y
# x⋅y⋅ℯ   + ℯ

# 不定積分: ∫ cos(x) dx
indef_integral = sympy.integrate(sympy.cos(x), x)
print("\n不定積分:")
sympy.pprint(indef_integral)
# 出力:
# 不定積分:
# sin(x)

# 定積分: ∫[0, pi] sin(x) dx
def_integral = sympy.integrate(sympy.sin(x), (x, 0, sympy.pi))
print("\n定積分:")
sympy.pprint(def_integral)
# 出力:
# 定積分:
# 2

# テイラー級数展開: e^x を x=0 の周りで4次まで展開
series_expr = sympy.exp(x).series(x, 0, 4)
print("\nテイラー級数展開:")
sympy.pprint(series_expr)
# 出力:
# テイラー級数展開:
#     2    3
#    x    x      ⎛ 4⎞
# 1 + x + ── + ── + O⎝x ⎠
#         2    6

方程式の求解

sympy.solve() 関数は、代数方程式、連立方程式、簡単な微分方程式などを解くための強力なツールです。

• 代数方程式: sympy.solve(方程式または式, 解く変数) で解を求めます。式を与えた場合は、式 = 0 として解かれます。方程式は sympy.Eq(左辺, 右辺) で表現できます。
• 連立方程式: 方程式のリストと解く変数のリストを渡します。
• 微分方程式: sympy.dsolve(微分方程式, 解く関数) を使用します。関数は sympy.Function で定義します。

sympy.init_printing(use_unicode=True)
x, y, a, b, c = sympy.symbols('x y a b c')

# 2次方程式: a*x**2 + b*x + c = 0 を x について解く
eq1 = sympy.Eq(a*x**2 + b*x + c, 0)
sol1 = sympy.solve(eq1, x)
print("2次方程式の解:")
sympy.pprint(sol1)
# 出力:
# 2次方程式の解:
# ⎡        ___________   ⎛       ___________⎞⎤
# ⎢       ╱           2  ⎜      ╱           2 ⎟⎥
# ⎢-b + ╲╱  -4⋅a⋅c + b   ⎜-b - ╲╱  -4⋅a⋅c + b  ⎟⎥
# ⎢─────────────────────, ⎜─────────────────────⎟⎥
# ⎣         2⋅a          ⎝         2⋅a         ⎠⎦


# 連立方程式:
# x + y = 2
# x - y = 0
eq2_1 = sympy.Eq(x + y, 2)
eq2_2 = sympy.Eq(x - y, 0)
sol2 = sympy.solve([eq2_1, eq2_2], (x, y))
print("\n連立方程式の解:")
sympy.pprint(sol2)
# 出力:
# 連立方程式の解:
# {x: 1, y: 1}

# 微分方程式: f'(x) = f(x)
f = sympy.Function('f')
eq3 = sympy.Eq(f(x).diff(x), f(x))
sol3 = sympy.dsolve(eq3, f(x))
print("\n微分方程式の解:")
sympy.pprint(sol3)
# 出力:
# 微分方程式の解:
#           x
# f(x) = C₁⋅ℯ

線形代数（行列）

sympy.Matrix クラスを使って行列を定義し、様々な操作を行うことができます。

• 行列の作成: リストのリストとして要素を指定します。
• 基本演算: 加算、減算、乗算（*）、べき乗（**）が可能です。
• 行列式: .det() メソッドで計算します。
• 逆行列: .inv() メソッドまたは **(-1) で計算します。
• 転置行列: .T 属性で取得します。
• 固有値・固有ベクトル: .eigenvals(), .eigenvects() メソッドで計算します。

sympy.init_printing(use_unicode=True)

# 行列の作成
M = sympy.Matrix([[1, 2], [3, 4]])
N = sympy.Matrix([[0, 1], [1, 0]])
vec = sympy.Matrix([x, y])

print("行列 M:")
sympy.pprint(M)
# 出力:
# 行列 M:
# ⎡1  2⎤
# ⎢    ⎥
# ⎣3  4⎦

# 行列演算
print("\nM + N:")
sympy.pprint(M + N)
# 出力:
# M + N:
# ⎡1  3⎤
# ⎢    ⎥
# ⎣4  4⎦

print("\nM * N:")
sympy.pprint(M * N)
# 出力:
# M * N:
# ⎡2  1⎤
# ⎢    ⎥
# ⎣4  3⎦

print("\nM * vec:")
sympy.pprint(M * vec)
# 出力:
# M * vec:
# ⎡ x + 2⋅y ⎤
# ⎢         ⎥
# ⎣3⋅x + 4⋅y⎦

# 行列式
det_M = M.det()
print(f"\ndet(M): {det_M}") # 出力: det(M): -2

# 逆行列
inv_M = M.inv()
print("\nM の逆行列:")
sympy.pprint(inv_M)
# 出力:
# M の逆行列:
# ⎡-2   1 ⎤
# ⎢       ⎥
# ⎣3/2  -1/2⎦

# 固有値
eigenvalues = M.eigenvals() # 固有値をキー、その重複度を値とする辞書を返す
print("\nM の固有値:")
sympy.pprint(eigenvalues)
# 出力:
# M の固有値:
# ⎧5   √33       √33   5   ⎫
# ⎨─ + ───: 1, - ─── + ─: 1⎬
# ⎩2    2         2    2   ⎭

# 固有ベクトル
eigenvectors = M.eigenvects() # (固有値, 重複度, [固有ベクトル]) のリストを返す
print("\nM の固有ベクトル:")
for val, mult, vec_list in eigenvectors:
    print("固有値:", val)
    print("固有ベクトル:", vec_list)
# 出力例（表示は環境によって多少異なる場合があります）:
# M の固有ベクトル:
# 固有値: 5/2 - sqrt(33)/2
# 固有ベクトル: [Matrix([
# [-2/3 - (-sqrt(33)/2 + 5/2)/3],
# [                             1]])]
# 固有値: 5/2 + sqrt(33)/2
# 固有ベクトル: [Matrix([
# [-2/3 - (sqrt(33)/2 + 5/2)/3],
# [                            1]])]

プロット

SymPyには基本的なプロット機能も含まれています（sympy.plottingモジュール）。簡単な関数のグラフを描画するのに便利です。より高度なプロットにはMatplotlibなどの専門ライブラリと連携させることが多いです。

from sympy.plotting import plot

p1 = plot(x**2, (x, -5, 5), title='y = x^2', show=False)
p2 = plot(sympy.sin(x), (x, -2*sympy.pi, 2*sympy.pi), line_color='red', show=False)

# 複数のプロットを重ねて表示
p1.extend(p2)
p1.show()

# パラメトリックプロット (媒介変数表示)
# sympy.plotting.plot_parametric((x座標の式, y座標の式), (媒介変数, 開始値, 終了値))
plot_parametric = sympy.plotting.plot_parametric((sympy.cos(x), sympy.sin(x)), (x, 0, 2*sympy.pi), title='円', aspect_ratio=(1,1))

# 陰関数プロット
# sympy.plotting.plot_implicit(方程式)
# plot_implicit = sympy.plotting.plot_implicit(sympy.Eq(x**2 + y**2, 5))

注意: SymPyのプロット機能を利用するには、バックエンドとなるライブラリ（通常はMatplotlib）がインストールされている必要があります。 show=False を指定するとプロットオブジェクトが返され、後で表示したり、他のプロットと組み合わせたりできます。

応用例と高度な機能

物理学や工学での利用

SymPyは、複雑な数式を扱う必要がある物理学や工学の分野で広く利用されています。

• 力学: 運動方程式の導出、ラグランジアンやハミルトニアンの計算など。
• 電磁気学: マクスウェル方程式の解析、ベクトル解析など。
• 量子力学: シュレディンガー方程式の解析、演算子の交換関係の計算など。
• 制御工学: 伝達関数の操作、ラプラス変換、システムの安定性解析など。
• 一般相対性理論: テンソル計算（EinsteinPyなど、SymPyを利用する外部ライブラリと連携）。

sympy.physics モジュールには、ベクトル演算、古典力学、量子力学などのためのサブモジュールが含まれています。

from sympy.physics.vector import ReferenceFrame, dynamicsymbols
from sympy import symbols

# 基準座標系の定義
N = ReferenceFrame('N')

# 時間変化するシンボル（動的シンボル）
theta = dynamicsymbols('theta')
omega = dynamicsymbols('omega') # 角速度
alpha = dynamicsymbols('alpha') # 角加速度

# 関係性の定義 (例: omega = d(theta)/dt)
# dynamicsymbols は時間 t の関数として扱われる
omega_def = theta.diff(symbols('t'))
print(f"ω = {omega_def}") # 出力: ω = Derivative(theta(t), t)

# ベクトル表現
# 例えば、x軸方向の単位ベクトルは N.x
Position = x*N.x + y*N.y

# 時間微分
Velocity = Position.dt(N) # N 基準系での時間微分
print(f"\n速度ベクトル: {Velocity}")
# 出力: 速度ベクトル: Derivative(x(t), t)*N.x + Derivative(y(t), t)*N.y

コード生成

SymPyで導出した数式を、NumPy、SciPy、C、Fortran、JavaScriptなどの他の言語やライブラリで使用できるコードに変換する機能があります（sympy.utilities.codegen モジュールや lambdify 関数）。これにより、記号計算で得られた複雑な数式を、高速な数値計算コードに組み込むことが容易になります。

from sympy import lambdify
import numpy as np

x, y = sympy.symbols('x y')
expr = sympy.sin(x) + sympy.cos(y)**2

# NumPy関数に変換
# lambdify(引数のリスト, 式)
func_numpy = lambdify((x, y), expr, 'numpy')

# NumPy配列を使って評価
x_vals = np.linspace(0, np.pi, 5)
y_vals = np.linspace(0, np.pi/2, 5)
result_numpy = func_numpy(x_vals, y_vals)

print("NumPy配列での評価結果:")
print(result_numpy)

# Cコードの生成例
from sympy.utilities.codegen import codegen
code_gen_result = codegen(('my_func', expr), 'C', 'my_module')
# code_gen_result には C のヘッダファイル (.h) とソースファイル (.c) の内容が含まれる
# print(code_gen_result[0][1]) # ヘッダファイルの内容
# print(code_gen_result[1][1]) # ソースファイルの内容

lambdify は、SymPyの式を効率的な数値計算用の関数に変換するのに特に便利です。Matplotlibでのプロットや、SciPyの最適化・積分ルーチンとの連携などで活躍します。

拡張性

SymPyは純粋なPythonで書かれているため、Pythonの豊富なエコシステムとシームレスに連携できます。また、ユーザーが独自の関数やオブジェクトを定義してSymPyの機能を拡張することも比較的容易です。多くの科学技術計算ライブラリ（例: PyDy (多体動力学), EinsteinPy (一般相対性理論), ChemPy (化学) など）がSymPyを基盤として利用、または連携して開発されています。

まとめ

SymPyは、Pythonで記号計算を行うための非常に強力で柔軟なライブラリです。基本的な代数演算から微積分、線形代数、方程式求解、さらにはコード生成まで、幅広い数学的操作を記号的に実行できます。

オープンソースであり、Pythonのエコシステムとの親和性も高いため、研究、教育、開発など様々な場面で活用されています。数式を手計算する手間を省き、より複雑な問題の解析に集中することを可能にしてくれます。

この記事ではSymPyの基本的な機能といくつかの応用例を紹介しましたが、SymPyにはここで紹介しきれなかった機能がまだまだたくさんあります。ぜひ公式ドキュメントなどを参照して、その可能性を探求してみてください。Let’s enjoy symbolic computation with SymPy!