【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]
下記は誤植と思われますので、ご確認ください。
「方程式・連立方程式を解く」の下の文で、解はタプル型 → 解はリスト型
修正しました。
ありがとうございます。