【SymPy】代数方程式

【SymPy】方程式・連立方程式を解く

SymPy の sympy.solve() を使って 代数方程式 $f(x)=0$ を解くことができます。

sympy.solve(f(x), x)

解はリスト型で得られます。例として $x^2+1=0$ を解いてみます。

# SYMPY_EQUATION_SOLVE_EQUATION

# In[1]

import sympy

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

# x**2 + 1 = 0 をxについて解く
s = sympy.solve(x**2 + 1, x)

print(s)
# [-I, I]

I は虚数単位を表す記号です。
すなわち $x^2+1=0$ の解は $x=i,\ -i$ です。
次は 3 次方程式を解いてみましょう。

# In[2]

# x**3-x**2-5*x+2=0をxについて解く
s = sympy.solve(x**3 - x**2 -5*x + 2, x)

print(s)
# [-2, -sqrt(5)/2 + 3/2, sqrt(5)/2 + 3/2]

連立方程式
 \[\begin{align*}&f(x)=0\\[6pt]&g(x)=0\\[6pt]&\cdots\cdots\cdots\end{align*}\]
の解は次の構文で取得することができます。

sympy.solve([f(x), g(x), ...], [x, y, ...])

# In[3]

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

# f(x)とg(x)を定義
fx = x + 2*y + 1
gx = 2*x + 5*y -3

# 連立方程式 f(x)=0, g(x)=0 の解
s = sympy.solve([fx, gx], [x, y])

print(s)
# {x: -11, y: 5}

円分方程式（円周等分方程式）

sympy.solve() で円分方程式 $x^n-1=0$ を解いてみましょう。円分方程式の解は複素数の範囲で必ず $n$ 個存在する（重解がない）ことが知られています。ここでは $n=5$ の場合、すなわち $x^5-1=0$ の解を求めてみます。ただし、得られる解の表式が複雑なので、MathJax 機能が組込まれた Jupyter Notebook の使用を前提とします。

# SYMPY_CYCLOTOMIC_EQUATION

# In[1]

import sympy

# 実行結果をLateXで表示
sympy.init_printing()

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

# f(x)を定義
fx = x**5 - 1

# 円分方程式 f(x)=0 の解
sympy.solve(fx, x)

実行結果は次のように LaTeX 形式で表示されます。
\[\left [ 1, \quad – \frac{1}{4} + \frac{\sqrt{5}}{4} – i \sqrt{\frac{\sqrt{5}}{8} + \frac{5}{8}}, \quad – \frac{1}{4} + \frac{\sqrt{5}}{4} + i \sqrt{\frac{\sqrt{5}}{8} + \frac{5}{8}},\\- \frac{\sqrt{5}}{4} – \frac{1}{4} – i \sqrt{- \frac{\sqrt{5}}{8} + \frac{5}{8}}, \quad – \frac{\sqrt{5}}{4} – \frac{1}{4} + i \sqrt{- \frac{\sqrt{5}}{8} + \frac{5}{8}}\right]\]

Eqオブジェクトによる方程式の定義

SymPyでは sympy.Eq() を使って方程式を定義することができます。

sympy.Eq(左辺, 右辺)

sympy.Eq() によって生成されたオブジェクトは sympy.core.relational.Equalityクラスに属しますが、名前が長いので以降は略して Eqオブジェクトとよぶことにします。sympy.Eq に渡した左辺と右辺の数式は lhs と rhs というデータ属性（インスタンス変数）に格納されるので、Eq.lhs, Eq.rhs によって左辺と右辺を取り出すことができます。sympy.solve() の引数に Eqオブジェクトと変数 var を渡すと lhs == rhs の解を返してきます。複雑な方程式を扱う場合は、sympy.solve() に数式を直接書きこむよりも、Eqオブジェクトを渡したほうがコードの可読性が高くなります。

# SYMPY_EQUATION_OBJECT

# In[1]

import sympy

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

# x**2 = 1
eq = sympy.Eq(x**2, 2*x + 1)

# 解を取得
s = sympy.solve(eq, x)

print("方程式の左辺 {}".format(eq.lhs))
print("方程式の右辺 {}".format(eq.rhs))
print("方程式の解 :", s)

# 方程式の左辺 x**2
# 方程式の右辺 2*x + 1
# 方程式の解 [1 + sqrt(2), -sqrt(2) + 1]