【SymPy】二次方程式
未知数の2乗が含まれる方程式を 二次方程式 とよびます。未知数が１個の場合、一般的には未知数を $x$、各項の係数を実数 $a,\ b,\ c$ として、
$$ax^2+bx+c=0\quad (a\neq 0)\tag{1}$$のような形で記述します。$b,\ c$ は $0$ でも構いませんが、$a$ が $0$ だと
$$bx+c=0\tag{2}$$のように一次以下の方程式になってしまうので、$a\neq 0$ は必須条件です。係数 $a,\ b,\ c$ の値によらず、二次方程式 (1) は複素数の範囲で必ず解をもつことが知られています。
それでは、簡単な二次方程式を解いてみましょう。(1) で $a=1,\ b=0, c=-9$ とします：
$$x^2-9=0\tag{3}$$左辺の $9$ を右辺に移項します。
$$x^2=9\tag{4}$$$x$ は $2$ 乗して $9$ になる数です。$x=3$ がすぐに頭に浮かぶと思いますが、符号を逆にした $-3$ も $2$ 乗すると $9$ になるので、(3) の解です。
二次方程式の解が複素数になることがあります。たとえば、(1) で $a=1,\ b=0, c=4$ とします：
$$x^2+4=0\tag{5}$$左辺の $4$ を右辺に移項すると
$$x^2=-4\tag{6}$$$2$ 乗して $-4$ なる数 $x$ は実数の範囲では存在しません。しかし、複素数の $2i$ と $-2i$ は (6) を満たします。すなわち $x=\pm 2i$ が (5) の解です。このような解を虚数解といいます。
SymPy を使って二次方程式 (5) を解いてみましょう。SymPy で代数方程式を解くときは、sympy.Eq で 二次方程式を定義（Equalityオブジェクトを作成）して、solve 関数に渡して解を得ます。
# In[1]
import sympy as sym
# 記号xを定義
x = sym.Symbol('x')
# 二次方程式x**2+2*x-3=0を定義
equation = sym.Eq(x**2 + 4, 0)
# 二次方程式をxについて解く
solutions = sym.solve(equation, x)
print(solutions)
# [-2*I, 2*I]
I は SymPy で虚数単位を表す記号です。もう少し複雑な二次方程式を考えてみましょう。
$$x^2+2x-3=0\tag{7}$$二次方程式は、左辺を因数分解できれば、比較的簡単に解を求められることが知られています。(7) の左辺を因数分解すると、
$$(x-1)(x+3)=0\tag{8}$$となります。因数 $x-1$ と $x+3$ のどちらがが $0$ ならば、因数の積 $(x-1)(x+3)$ は $0$ になるので、解は $x=1$ または $x=-3$ となります。念のために、SymPy で検算しておきましょう。
# In[2] import sympy as sym # 二次方程式x**2+2*x-3=0を定義 equation = sym.Eq(x**2 + 2*x -3, 0) # 二次方程式をxについて解く solutions = sym.solve(equation, x) print(solutions) # [-3, 1]
とはいえ、左辺が必ずしも簡単に因数分解できるとは限りません。むしろ、できないケースが大半です。たとえば、
$$11x^2-23x+8=0\tag{9}$$のような方程式の解は分数と根号を含むため、実質上、手計算で因数分解することは不可能です。このような場合、以下で説明する解の公式を用います。
二次方程式の解の公式
先にも述べたように、二次方程式
$$ax^2+bx+c=0\quad (a\neq 0)\tag{10}$$は複素数の範囲で必ず解をもちます。では、係数がすべて実数の時、(10) の一般解がどのような形をしているのか、SymPy で求めてみましょう。結果を見やすい形にするために、init_printing 関数で、LaTex 形式の数式を出力するように設定しておきます。
# In[3]
sym.init_printing()
# 実数係数a,b,cを定義
a, b, c = sym.symbols('a:c', real=True)
# 二次方程式f(x)=0のf(x)を定義
equation = sym.Eq(a*x**2 + b*x + c, 0)
# 二次方程式f(x)=0を解く
solutions = sym.solve(equation, x)
# 解を表示する
display(solutions)
これが、いわゆる 解の公式 とよばれるものです。根号の前にある符号を±でまとめると、
$$\frac{- b \pm \sqrt{b^{2}-4ac}}{2 a}\tag{11}$$となります。数学の教科書に必ず載っていますが、ここでも一応、手計算で (10) を導いてみましょう（面倒であれば、ここの部分は読み飛ばしてください）。(1) の左辺を $p(x)$ とおきます。
$$p(x)=ax^2+bx+c\tag{12}$$最初の２つの項を $a$ で括ります。
$$p(x)=a\left(x^2+\frac{b}{a}x\right)+c\tag{13}$$ここで、平方完成という手法を用います。括弧で括られた部分に着目して、無理に平方の形を作って
$$p(x)=a\left(x+\frac{b}{2a}\right)^2-k+c\tag{14}$$とおきます。ここで、$k$ は $\displaystyle a\left(x+\frac{b}{2a}\right)^2$ と $\displaystyle a\left(x^2+\frac{b}{a}\right)$ の差分です。実際、式 (13) 第１項目を展開してみると
$$a\left(x+\frac{b}{2a}\right)^2=a\left(x^2+\frac{b}{a}x+\frac{b^2}{4a^2}\right)\tag{15}$$
なので、$\displaystyle\frac{b^2}{4a}$ だけ余分な項が増えていて、(12) と一致させるためには、これを差し引かなければなりません。すなわち、$\displaystyle k=\frac{b^2}{4a}$ となります。
$$p(x)=a\left(x+\frac{b}{2a}\right)^2-\frac{b^2}{4a}+c=a\left(x+\frac{b}{2a}\right)^2-\frac{b^2-4ac}{4a}\tag{16}$$したがって、二次方程式 $p(x)=0$ は
$$a\left(x+\frac{b}{2a}\right)^2=\frac{b^2-4ac}{4a}\tag{17}$$となるので、両辺を $a$ で割って
$$\left(x+\frac{b}{2a}\right)^2=\frac{b^2-4ac}{4a^2}\tag{18}$$さらに両辺の平方根をとります。
$$x+\frac{b}{2a}=\pm\frac{\sqrt{b^2-4ac}}{2a}\tag{19}$$したがって、求める解 $x$ は
$$x=\frac{-b\pm\sqrt{b^2-4ac}}{2a}\tag{20}$$となります。
さきほどの二次方程式
$$11x^2-23x+8=0\tag{21}$$を解の公式を使って解いてみましょう。解の公式 (20) に $a=11,\ b=-23, c=8$ を代入して
$$x=\frac{23\pm\sqrt{23^2-4\times 11\times 8}}{2\times 11}=\frac{23\pm\sqrt{177}}{22}\tag{22}$$となります。SymPy で検算してみましょう。
# In[4] # 二次方程式11*x**2-23*x+8=0を定義 equation = sym.Eq(11*x**2 - 23*x + 8, 0) # 二次方程式をxについて解く solutions = sym.solve(equation, x) display(solutions) sym.latex(solutions)
二次方程式の判別式
解の公式
$$x=\frac{-b\pm\sqrt{b^2-4ac}}{2a}$$をよく眺めてみると、根号の中の $b^2-4ac$ の値によって、解が実数なのか複素数なのか判別できることがわかります。$D=b^2-4ac$ とおくと、
[1] $D$ が正のとき、二次方程式は2個の実数解をもちます。
[2] $D=0$ のときは $\displaystyle x=-\frac{b}{2a}\tag{21}$ となるので、ただ１個の実数解をもつことがわかります。ちなみに、このような解を重解とよびます。
[3] $D$ が負の時は２個の虚数解をもちます。
$D=b^2-4ac$ を二次方程式の 判別式 とよびます。$D$ は英語の discriminant の頭文字です。SymPy には判別式の値を計算する discriminant 関数が用意されています。さきほどの二次方程式 (21)
$$11x^2-23x+8=0$$の判別式を計算してみましょう。
# In[5] # 二次方程式の左辺の式を使って判別式の値を計算 D = sym.discriminant(equation.lhs) print(D) # 177
解 $x$ を実数に制限した場合、$D\lt 0$ は「解なし」です。SymPy では、記号を定義する際に real=True を設定することで記号を実数に限定できます。この条件の下で、たとえば
$$x^2+2x+10=0\tag{23}$$を SymPy に解かせてみると空白のタプルを返します。
# In[6]
x = sym.Symbol('x', real=True)
# 二次方程式f(x)=0のf(x)を定義
equation = sym.Eq(x**2 + 2*x + 10, 0)
# 二次方程式f(x)=0を解く
solutions = sym.solve(equation, x)
display(solutions)
#[]
二次関数のグラフを使って解の個数を調べる
二次方程式 $ax^2+bx+c=0$ の解は二次関数 $f(x)=ax^2+bx+c$ と $g(x)=0$ の交点とみなすことができます。以下のコードで、係数を変えた３つの二次関数のグラフを並べて表示します。
# In[7]
import numpy as np
import matplotlib.pyplot as plt
# xの離散値リストを作成
x = np.linspace(-5, 5, 100)
# 3行1列のサブプロットを用意する
fig, axes = plt.subplots(nrows=3, ncols=1, figsize=(4, 10))
# 二次関数の係数の組み合わせ
coefficients = [(2, 5, 1), (1, 2, 1), (1, 2, 5)]
for i, ax in enumerate(axes.flat):
# グラフタイトルに係数と判別式の値を表示
ax.set_title(f"a={a}, b={b}, c={c}, D={d}")
# y=0の位置にx軸を太線で引く
ax.axhline(y=0, color='black', linewidth=2, zorder=1)
# 軸範囲と軸ラベル、グリッド線の表示
ax.set_xlim(-5, 5)
ax.set_ylim(-4, 8)
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.grid(True)
# 二次関数の値と判別式の値を計算
a, b, c = coefficients[i]
y = a * x**2 + b * x + c
d = b**2 - 4*a*c
# 二次関数をプロット
ax.plot(x, y, zorder=0)
# x軸との交点を表示
if d == 0:
solution = -b/(2*a)
ax.scatter(solution, 0, c='red', zorder=2)
elif d > 0:
solution = [(-b+d**0.5)/(2*a), (-b-d**0.5)/(2*a)]
ax.scatter(solution, [0, 0], c='red', zorder=2)
plt.tight_layout()
plt.show()
上のグラフからも明らかなように、$D\gt 0$ のときは２個の交点、$D=0$ のときは１個の交点（接点）をもち、そして $D\lt 0$ のときは交点をもちません。
解と係数の関係
二次方程式の解を
$$\alpha=\frac{- b+\sqrt{b^2-4ac}}{2 a},\quad \beta=\frac{- b-\sqrt{b^2-4ac}}{2 a}\tag{24}$$とおいて、和と積を作ってみると、
$$\alpha+\beta=-\frac{b}{a},\quad \alpha\beta=\frac{c}{a}\tag{25}$$となります。つまり、二次方程式を解かなくても、$\alpha+\beta$ と $\alpha\beta$ は計算できます。二次方程式
$$3x^2+11x+5=0\tag{26}$$について、$\displaystyle\frac{1}{\alpha}+\frac{1}{\beta}$ を求めてみましょう。解と係数の関係により、
$$\alpha+\beta=-\frac{11}{3},\quad\alpha\beta=\frac{5}{3}\tag{27}$$なので、
$$\frac{1}{\alpha}+\frac{1}{\beta}=\frac{\alpha+\beta}{\alpha\beta}=-\frac{11}{5}\tag{28}$$となります。SymPy で検算してみます。
# In[8]
x = sym.Symbol('x')
# 二次方程式f(x)=0のf(x)を定義
equation = sym.Eq(3*x**2 + 11*x + 5, 0)
# 二次方程式f(x)=0を解く
alpha, beta = sym.solve(equation, x)
# 1/alpha+1/beta
print(sym.simplify(1/alpha + 1/beta))
# -11/5
ランダム係数二次方程式
SymPy を活用したサンプルコードを作っておきました。後半に紹介する Rquadratic() は二次方程式などを自動的に作成するクラスなので、数学の試験問題作成などに役立ててください。
二次方程式クラス作成の準備として、有理係数二次多項式クラス (Quadratic) を定義しておきます。Quadratic は二次多項式 $ax^2+bx+c$ をデータ属性にもつインスタンスを生成します。各種のメソッドによって x に値を代入したり、頂点の座標や二次方程式 $ax^2+bx+c=0$ の解を求めることができます。
# PYTHON_QUADRATIC_EQUATION
# In[1]
import sympy
# 有利係数2次多項式クラス
class Quadratic:
def __init__(self, a, b, c):
# 記号xを定義
sympy.var('x')
# 係数a,b,cを有理数型に変換
self.a = sympy.Rational(a)
self.b = sympy.Rational(b)
self.c = sympy.Rational(c)
# 2次多項式の表式
self.exps = self.a*x**2 + self.b*x + self.c
# 2次多項式の値を返すメソッド
def value(self, t):
return self.exps.subs(x, t)
# 頂点の座標を返すメソッド
def vertex(self):
xv = -self.b / (2*self.a)
yv = (-self.b**2 + 4*self.a*self.c)/(4*self.a)
return(xv, yv)
# 2次方程式と解を返すメソッド
def equation(self):
eq = sympy.Eq(self.exps, 0)
return (eq, sympy.solve(eq))
Quadratic(a, b, c) は二次多項式オブジェクト $ax^2+bx+c$ を生成します。渡された係数 a, b, c は sympy.Rational() によってすべて有理数型に変換されます（たとえば 1.5 は 3/2 に変換されます）。
Quadratic.value(p) は二次多項式オブジェクトの記号 x に p の値を代入した結果を返します。
Quadratic.vertex() は二次関数の頂点の座標をタプルで返します。
Quadratic.equation() は二次方程式 $ax^2+bx+c=0$ の表式と解をタプルで返します。解はタプルの中にリスト形式で収められています。
以下に Quadraticクラスの使用例を載せておきます。
# In[2] # Quadraticクラスのインスタンスを生成 y = Quadratic(1, 2, 3) # 2次多項式x**2+2*x+3 print(y.exps) # 2次関数y=x**2+2*x+3 の頂点の座標 print(y.vertex()) # 2次方程式x**2+2*x+3=0 と解の表示 print(y.equation()) # x**2 + 2*x + 3 # (-1, 2) # (Eq(x**2 + 2*x + 3, 0), [-1 - sqrt(2)*I, -1 + sqrt(2)*I])
Jupyter Notebook で実行する場合は、以下のコードを入力すると結果が LaTeX 数式で表示されます (display() は Jupyter Notebook の関数です)。
sympy.init_printing() # Jupyter Notebook LaTeX表示バージョン # Quadraticクラスのインスタンスを生成 y = Quadratic(1, 2, 3) # 2次多項式 x**2+2*x+3 display(y.exps) # x**2 + 2*x + 3 = 0 と解の表示 display(y.equation())
&\left ( x^{2} + 2 x + 3 = 0, \quad \left [ -1 – \sqrt{2} i, \quad -1 + \sqrt{2} i\right ]\right )\end{align*}\]
それでは、ランダム係数二次方程式を生成する Rquadratic クラスを定義してみます。Rquadratic は Quadratic のサブクラスです。その機能のほとんどを Quadratic から継承していますが、係数 $a,\ b,\ c$ はランダムに決定される整数です。
# In[3]
import random
# ランダム係数2次方程式クラス
class Rquadratic(Quadratic):
# 初期化メソッドのオーバーライド
def __init__(self, abc_range):
# 記号xを定義
sympy.var('x')
# 2次多項式の係数a,b,cをランダムに決定
self.a = random.randint(abc_range[0], abc_range[1])
self.b = random.randint(abc_range[0], abc_range[1])
self.c = random.randint(abc_range[0], abc_range[1])
# 2次多項式の表式
self.exps = self.a*x**2 + self.b*x + self.c
Rquadratic(abc_range) はランダムな整数係数をもつ二次多項式オブジェクトを生成します。引数 abc_range は各係数の下限値と上限値です。
Rquadratic.equation() は二次方程式 $ax^2+bx+c=0$ の表式と解をタプルで返します。解はタプルの中にリスト形式で収められます。
以下に Rquadraticクラスの使用例を載せておきます。
# In[4]
random.seed(0)
# Rquadraticクラスのインスタンスを生成
# 2次多項式の係数a,b,cの範囲は[-5,5]
y = Rquadratic((-5, 5))
# ランダム係数の2次多項式
print(f'2次多項式：{y.exps}')
eq, sol = y.equation()
# 2次方程式
print(f'2次方程式：{eq}')
# 2次方程式の解
print(f'2次方程式の解：{sol}')
# 2次多項式：x**2 + x - 5
# 2次方程式：Eq(x**2 + x - 5, 0)
# 2次方程式の解：[-1/2 + sqrt(21)/2, -sqrt(21)/2 - 1/2]
コメント
下記は誤植と思われますので、ご確認ください。
「Quadraticクラスの定義」プログラムで、return [eq, sympy.solve(eq)] → return (eq, sympy.solve(eq))
修正しました。
ありがとうございます。m(_ _)m
【お知らせ】二次方程式を生成する Rquadratic クラスの定義文を見直してみると、無作為に定まる係数が正の値しかとれないことに気づいたので修正しました。修正された Rquadratic() のコンストラクタの引数 abc_range に、タプルやリストで (-3, 5) や [-3, 5] のように指定すると、係数 a, b, c が -5 ～ 3 の範囲で生成されます。また、記事全体を見直して細かな修正をしておきました。