sympy-symbolic-math

SymPy — 記号数学

概要

SymPy は、数値近似ではなく数学記号を使って厳密計算を行う、記号数学向けの Python ライブラリです。純粋な Python で、外部依存関係なしに、代数、微積分、方程式の求解、線形代数、物理学、コード生成を扱えます。

使用する場面

• 方程式を記号的に解く場合（代数方程式、連立方程式、微分方程式）
• 微積分の演算を行う場合（導関数、積分、極限、級数展開）
• 代数式を簡約・操作する場合
• 行列を記号的に扱う場合（固有値、行列式、分解）
• 記号式を高速な数値関数に変換する場合（lambdify → NumPy）
• 数式からコードを生成する場合（C、Fortran、LaTeX）
• 厳密な結果が必要な場合（例：1.414... ではなく sqrt(2)）
• 数値計算（配列演算、データに対する線形代数）には、numpy/scipy を使用します
• 統計モデリング（回帰、仮説検定）には、statsmodels を使用します

前提条件

pip install sympy
# Optional for numerical evaluation:
pip install numpy matplotlib

SymPy は純粋な Python で実装されているため、コンパイル済みの依存関係はなく、どこにでもインストールできます。

クイックスタート

from sympy import symbols, solve, diff, integrate, sqrt, pi

x = symbols('x')

# Solve equation
print(solve(x**2 - 5*x + 6, x))          # [2, 3]

# Derivative
print(diff(x**3 + 2*x, x))                # 3*x**2 + 2

# Integral
print(integrate(x**2, (x, 0, 1)))         # 1/3

# Exact arithmetic
print(sqrt(8))                              # 2*sqrt(2)
print(pi.evalf(30))                        # 3.14159265358979323846264338328

コア API

1. 記号と式

記号変数を作成し、式を操作します。

from sympy import symbols, Symbol, Rational, S, oo, pi, E, I
from sympy import simplify, expand, factor, collect, cancel, trigsimp

# Define symbols
x, y, z = symbols('x y z')

# With assumptions (improve simplification)
n = symbols('n', integer=True)
t = symbols('t', positive=True, real=True)
from sympy import sqrt
print(sqrt(t**2))   # t (not Abs(t), because t is positive)

# Exact fractions (avoid floats!)
expr = Rational(1, 3) * x + S(1)/7
print(expr)  # x/3 + 1/7

# Simplification
print(simplify(x**2 + 2*x + 1))          # (x + 1)**2
print(expand((x + 1)**3))                 # x**3 + 3*x**2 + 3*x + 1
print(factor(x**3 - x))                   # x*(x - 1)*(x + 1)
print(collect(x*y + x - 3 + 2*x**2 - z*x**2, x))  # x**2*(2 - z) + x*(y + 1) - 3

2. 微積分

導関数、積分、極限、級数を扱います。

from sympy import symbols, diff, integrate, limit, series, oo, sin, cos, exp, log

x = symbols('x')

# Derivatives
print(diff(sin(x**2), x))                 # 2*x*cos(x**2)
print(diff(x**4, x, 3))                   # 24*x (third derivative)

# Partial derivatives
x, y = symbols('x y')
f = x**2 * y**3
print(diff(f, x, y))                      # 6*x*y**2

# Integrals
x = symbols('x')
print(integrate(x**2, x))                 # x**3/3 (indefinite)
print(integrate(exp(-x**2), (x, -oo, oo)))  # sqrt(pi) (Gaussian)
print(integrate(x * exp(-x), (x, 0, oo)))  # 1

# Limits
print(limit(sin(x)/x, x, 0))             # 1
print(limit((1 + 1/x)**x, x, oo))        # E

# Taylor series
print(series(exp(x), x, 0, 5))           # 1 + x + x**2/2 + x**3/6 + x**4/24 + O(x**5)

3. 方程式を解く

代数方程式、超越方程式、微分方程式を解きます。

from sympy import symbols, solve, solveset, Eq, S, linsolve, nonlinsolve, Function, dsolve

x, y = symbols('x y')

# Single equation
print(solve(x**2 - 4, x))                 # [-2, 2]
print(solveset(x**2 - 4, x, S.Reals))     # {-2, 2}

# System of linear equations
print(linsolve([x + y - 5, 2*x - y - 1], x, y))  # {(2, 3)}

# System of nonlinear equations
print(nonlinsolve([x**2 + y - 4, x + y**2 - 4], x, y))

# Differential equation: y'' + y = 0
f = Function('f')
ode = f(x).diff(x, 2) + f(x)
print(dsolve(ode, f(x)))                  # Eq(f(x), C1*sin(x) + C2*cos(x))

# With initial conditions
from sympy import Derivative
ics = {f(0): 1, f(x).diff(x).subs(x, 0): 0}
print(dsolve(ode, f(x), ics=ics))         # Eq(f(x), cos(x))

4. 行列と線形代数

記号的な行列演算を行います。

from sympy import Matrix, eye, zeros, ones, diag, symbols

# Create matrices
M = Matrix([[1, 2], [3, 4]])
print(f"Det: {M.det()}")                   # -2
print(f"Inverse:\n{M**-1}")

# Symbolic matrices
a, b = symbols('a b')
M = Matrix([[a, b], [b, a]])
print(f"Eigenvalues: {M.eigenvals()}")     # {a - b: 1, a + b: 1}

# Eigenvectors and diagonalization
eigendata = M.eigenvects()
# [(eigenval, multiplicity, [eigenvectors]), ...]
P, D = M.diagonalize()
print(f"M = P*D*P^-1")

# Solve linear system Ax = b
A = Matrix([[1, 2], [3, 4]])
b = Matrix([5, 6])
x = A.solve(b)
print(f"Solution: {x.T}")

# Matrix calculus
t = symbols('t')
M_t = Matrix([[t, t**2], [1, t]])
print(f"dM/dt:\n{M_t.diff(t)}")

5. コード生成

記号式を高速な数値関数またはコンパイル済みコードに変換します。

import numpy as np
from sympy import symbols, lambdify, sin, exp, ccode, fcode, latex

x, y = symbols('x y')
expr = sin(x) * exp(-x**2 / 2)

# lambdify: symbolic → fast NumPy function
f = lambdify(x, expr, 'numpy')
x_vals = np.linspace(-5, 5, 1000)
y_vals = f(x_vals)
print(f"Shape: {y_vals.shape}, Max: {y_vals.max():.4f}")

# Multi-variable lambdify
expr2 = x**2 + y**2
f2 = lambdify((x, y), expr2, 'numpy')
print(f"f(3, 4) = {f2(3, 4)}")            # 25

# C code generation
print(ccode(expr))                          # sin(x)*exp(-1.0/2.0*pow(x, 2))

# Fortran code generation
print(fcode(expr))

# LaTeX output
print(latex(expr))                          # \sin{\left(x \right)} e^{- \frac{x^{2}}{2}}

6. 物理学モジュール

古典力学、ベクトル解析、単位を扱います。

from sympy import symbols, cos, sin, Function
from sympy.physics.mechanics import dynamicsymbols, LagrangesMethod, Particle, Point, ReferenceFrame
from sympy.physics.vector import dot, cross

# Vector analysis
N = ReferenceFrame('N')
v1 = 3*N.x + 4*N.y + 0*N.z
v2 = 1*N.x + 0*N.y + 2*N.z
print(f"Dot: {dot(v1, v2)}")               # 3
print(f"Cross: {cross(v1, v2)}")            # 8*N.x - 6*N.y - 4*N.z

# Simple pendulum via Lagrangian mechanics
q = dynamicsymbols('q')       # Generalized coordinate (angle)
m, g, l = symbols('m g l', positive=True)
T = Rational(1, 2) * m * (l * q.diff())**2          # Kinetic energy
V = m * g * l * (1 - cos(q))                         # Potential energy
L = T - V                                            # Lagrangian
print(f"Lagrangian: {L}")

重要な概念

厳密演算と数値演算

from sympy import Rational, S, sqrt, pi

# WRONG: introduces floating-point error
expr_bad = 0.5 * x           # Float 0.5, loses exactness

# CORRECT: exact symbolic arithmetic
expr_good = Rational(1, 2) * x    # Exact 1/2
expr_good = S(1)/2 * x            # Alternative exact syntax
expr_good = x / 2                  # Also exact

# Numerical evaluation when needed
print(sqrt(2).evalf())         # 1.41421356237310
print(pi.evalf(50))            # 50 digits of precision

ソルバーの選択ガイド

ソルバー	使用する場面	戻り値
solve(eq, x)	汎用、レガシー	解のリスト
solveset(eq, x, domain)	代数方程式（推奨）	集合（無限集合の場合あり）
linsolve(system, vars)	線形系	タプルの FiniteSet
nonlinsolve(system, vars)	非線形系	タプルの FiniteSet
dsolve(ode, f(x))	常微分方程式	等式（Eq）
nsolve(eq, x0)	数値的な根の探索	浮動小数点数による近似

よく使われる簡約関数

関数	処理内容	例
simplify()	一般的な簡約（低速で、あらゆる方法を試します）	sin(x)**2 + cos(x)**2 → 1
expand()	乗算を分配	(x+1)**2 → x**2+2*x+1
factor()	既約因子に因数分解	x**2-1 → (x-1)*(x+1)
collect()	変数ごとにまとめる	x の項をまとめます
cancel()	分数の共通因子を約分	(x**2-1)/(x-1) → x+1
trigsimp()	三角関数の式を簡約	三角関数では simplify より高速
powsimp()	べき乗・指数関数を簡約	x**a * x**b をまとめます

よくあるワークフロー

ワークフロー：記号式から数値計算へのパイプライン

from sympy import symbols, diff, integrate, lambdify, sin, cos
import numpy as np
import matplotlib.pyplot as plt

x = symbols('x')

# 1. Define expression symbolically
f_expr = sin(x) * cos(x)**2

# 2. Symbolic operations
f_prime = diff(f_expr, x)
F_expr = integrate(f_expr, x)
print(f"f(x) = {f_expr}")
print(f"f'(x) = {f_prime}")
print(f"F(x) = {F_expr}")

# 3. Convert to fast numerical functions
f_num = lambdify(x, f_expr, 'numpy')
f_prime_num = lambdify(x, f_prime, 'numpy')
F_num = lambdify(x, F_expr, 'numpy')

# 4. Evaluate and plot
x_vals = np.linspace(0, 2*np.pi, 500)
fig, axes = plt.subplots(1, 3, figsize=(12, 4))
axes[0].plot(x_vals, f_num(x_vals)); axes[0].set_title('f(x)')
axes[1].plot(x_vals, f_prime_num(x_vals)); axes[1].set_title("f'(x)")
axes[2].plot(x_vals, F_num(x_vals)); axes[2].set_title('F(x)')
plt.tight_layout()
plt.savefig('symbolic_pipeline.png', dpi=150)
print("Saved symbolic_pipeline.png")

ワークフロー：解いて検証する

from sympy import symbols, solve, simplify, Eq, sqrt

x = symbols('x')

# 1. Define equation
equation = x**3 - 6*x**2 + 11*x - 6

# 2. Solve symbolically
solutions = solve(equation, x)
print(f"Solutions: {solutions}")            # [1, 2, 3]

# 3. Verify each solution
for sol in solutions:
    result = simplify(equation.subs(x, sol))
    assert result == 0, f"Solution {sol} failed!"
    print(f"  x={sol}: f(x) = {result} ✓")

# 4. Factor the polynomial
from sympy import factor
print(f"Factored: {factor(equation)}")      # (x - 1)*(x - 2)*(x - 3)

ワークフロー：ODE 系の解析

1. Function と dsolve() を使って ODE を定義します
2. 記号的に解き、ics={} パラメーターで初期条件を適用します
3. lambdify() で解を数値関数に変換します
4. matplotlib で解の軌跡をプロットします

主要なパラメーター

パラメーター	関数	デフォルト	選択肢	効果
domain	solveset()	S.Complexes	S.Reals、S.Integers	解の領域を制限します
force	simplify()	False	True/False	積極的に簡約します
n	diff(expr, x, n)	1	1–∞	導関数の次数
精度	evalf(n)	15	1–1000+	数値精度の桁数
バックエンド	lambdify()	"math"	"numpy"、"scipy"、"mpmath"	評価に使用する数値バックエンド
rational	nsimplify()	True	True/False	厳密な有理数近似を求めます

ベストプラクティス

1. 分数には常に Rational() または S() を使用します：0.5 * x は浮動小数点数を導入し、厳密な計算に支障をきたします。Rational(1, 2) * x または S(1)/2 * x を使用してください。

2. シンボルに仮定を追加します：symbols('x', positive=True) を指定すると、sqrt(x**2) → x のような簡約が可能になります。仮定がない場合、SymPy は一般的な複素数のケースを扱う必要があります。

3. 数値評価には subs().evalf() ではなく lambdify を使用します：ループ内での subs/evalf は、lambdify の呼び出し 1 回と比べて 100～1000 倍遅くなります。

   # Slow: [expr.subs(x, v).evalf() for v in values]
   # Fast: f = lambdify(x, expr, 'numpy'); f(np.array(values))
   

4. アンチパターン — デフォルトで simplify() を使用する：simplify() は多くの方法を試すため、処理が遅くなります。目的の形式が分かっている場合は、専用の関数（factor、expand、trigsimp）を使用してください。

5. 代数方程式には solve より solveset を優先します：solveset は適切な数学的集合を返し、エッジケースにもより適切に対応します。solve はレガシーですが、一般的なケースでは引き続き便利です。

6. アンチパターン — 数値解で十分な場合に記号的に解く：閉形式の解がない方程式では、solve が失敗するのを待つのではなく、数値的な根を求めるために nsolve(eq, x0) を使用してください。

7. Jupyter では init_printing() を使用して出力を読みやすくします：from sympy import init_printing; init_printing() を実行すると、ノートブックで LaTeX 表示が有効になります。

よく使うレシピ

レシピ：LaTeX ドキュメントを生成する

from sympy import symbols, Integral, Eq, latex, sqrt, pi

x = symbols('x')
integral = Integral(x**2 * sqrt(1 - x**2), (x, 0, 1))
result = integral.doit()

print(f"$$ {latex(integral)} = {latex(result)} $$")
# $$ \int\limits_{0}^{1} x^{2} \sqrt{1 - x^{2}}\, dx = \frac{\pi}{16} $$

レシピ：パラメーター付き ODE の解

from sympy import symbols, Function, dsolve, Eq, exp, lambdify
import numpy as np

x = symbols('x')
k, A = symbols('k A', positive=True)
f = Function('f')

# Solve y' = -ky with y(0) = A
ode = Eq(f(x).diff(x), -k * f(x))
solution = dsolve(ode, f(x), ics={f(0): A})
print(f"Solution: {solution}")              # f(x) = A*exp(-k*x)

# Evaluate for specific parameters
f_num = lambdify((x, k, A), solution.rhs, 'numpy')
x_vals = np.linspace(0, 5, 100)
y_vals = f_num(x_vals, k=0.5, A=10)
print(f"y(5) = {y_vals[-1]:.4f}")

レシピ：記号行列の分解

from sympy import Matrix, symbols, pprint

a, b, c, d = symbols('a b c d')
M = Matrix([[a, b], [c, d]])

# Characteristic polynomial
lam = symbols('lambda')
char_poly = M.charpoly(lam)
print(f"Characteristic polynomial: {char_poly.as_expr()}")

# Eigenvalues (symbolic)
eigenvals = M.eigenvals()
print(f"Eigenvalues: {eigenvals}")

# Determinant and trace
print(f"det(M) = {M.det()}")               # a*d - b*c
print(f"tr(M) = {M.trace()}")              # a + d

トラブルシューティング

問題	原因	解決策
NameError: name 'x' is not defined	シンボルが作成されていない	使用する前に x = symbols('x') で定義する
予期しない浮動小数点数の結果	Rational(1,2) の代わりに 0.5 を使用している	厳密な分数には Rational() または S() を使用する
simplify() が非常に遅い	複雑な式に対してすべての方法を試している	専用の関数 factor()、expand()、trigsimp() を使用する
solve() が空のリストを返す	閉形式の解が存在しない	数値近似には nsolve(eq, x0) を使用する
sqrt(x**2) が x ではなく sqrt(x**2) を返す	x に仮定がない	x = symbols('x', positive=True) で定義する
lambdify の結果が間違っている	式に SymPy 固有の関数が含まれている	バックエンドを指定する：lambdify(x, expr, 'numpy') または 'scipy'
dsolve で NotImplementedError が発生する	ODE の型がサポートされていない	代わりに数値 ODE ソルバー（scipy の odeint）を試す

関連スキル

• matplotlib-scientific-plotting — lambdify による変換後に記号計算の結果をプロットする
• statsmodels-statistical-modeling — 統計的推論。厳密な代数ではなく p 値が必要な場合に使用する
• matlab-scientific-computing — 数値計算（記号計算ではない）向けの MATLAB の代替手段

参考資料

• SymPy documentation — 公式の API リファレンスとチュートリアル
• SymPy tutorial — 初心者向けの入門ガイド
• SymPy GitHub — ソースコード、例、課題
• Meurer et al. (2017)「SymPy: symbolic computing in Python」— PeerJ Computer Science