三次方程式を最も早く解ける解法はどれだ！？~Python編~
実は三次方程式にも解の公式が存在した！？
たまちゃん
三次方程式ってどうやって解くんだ？
有名なのは因数定理かな？解の公式もあるけど
えんちゃん
たまちゃん
え！三次方程式って解の公式あんのか！
おん。ちょいと複雑で高校では習わんかもな
えんちゃん
ということで今回は色んな手法を使って三次方程式の実数解を求めていきたいと思います。また、最終的にどの方法が一番早く解けるのかを調べていきます。ミソは実数解の取り出し方です。
目次
SympyとはPythonにおいて記号を用いて計算処理を行うことができるライブラリのことです。詳しい内容については前回の記事で書いているので参考にしてみてください。
例えば$$x-2=0$$を解きたい場合はPythonで以下のように書くことができます。
from sympy import *
x = symbols("x")
ans = solve(Eq(x-2, 0), x)
#[2]先ずx = symbols(“x”)の部分で変数xを定義します。そしてsolve(Eq(A,B),C)のA(左辺)、B(右辺)、C(解く対象)にそれぞれx-2、0、xを入れます。sympyの場合はこれで終わりです。
それではこれを例に実際に次のような三次方程式を解いてみましょう。$$x^3-3x^2+7x-5=0$$先ず先ほどの変数の定義などを行います。
def cubic_sympy(a,b,c,d):
x = symbols("x")
ans = solve(Eq(a * x ** 3 + b * x ** 2 + c * x + d, 0), x)]三次方程式を解くだけならこれで終わりです。しかし、もし実数解だけ欲しいという場合は虚数が含まれるかどうか判定する必要があります。実際に複数の答えの中から実数解だけを取り出すという処理を書いていきましょう。これには判別式を使ってもいいのですが、せっかくPythonを使っているのでここではライブラリを使った条件分岐を書きます。
x_list = []
for i in ans:
#虚数判定を行うためAdd型をcomplex型に直す
tmp = complex(i)
#実数解だけを抽出
if np.iscomplex(tmp) == False or np.iscomplex(i) and abs(tmp.imag) < 0.0000001:
x_list.append(tmp.real)
return x_list先ず3つの解を格納するリストx_listを作ります。そして先ほど求めた3つの解が格納されているansから1つずつ取り出して処理を行っていきます。後の虚数判定のために一度complex型に直してから条件分岐を書いていきます。具体的には解に虚数が含まれているかどうかを判定するnp.iscomplexを使用し、含まれていなかったらx_listに格納していきます。
また、虚数解が出ないはずなのに虚部が出ることがあります。これは誤差の表記として虚部が出ているだけで実際の解は実数なのでこれもx_listに格納します。そのために虚部の部分が限りなく小さい場合は格納するという風に書いています。ちなみにcomplex型の実部・虚部の取り出し方はそれぞれx.real、x.imagです。全体のコードは以下のようになります。
from sympy import *
import numpy as np
def cubic_sympy(a,b,c,d):
x = symbols("x")
ans = solve(Eq(a * x ** 3 + b * x ** 2 + c * x + d, 0), x)
x_list = []
for i in ans:
tmp = complex(i)
if np.iscomplex(tmp) == False or np.iscomplex(i) and abs(tmp.imag) < 0.0000001:
x_list.append(tmp.real)
return x_list
if __name__ == '__main__':
print(cubic_equation2(1,-2,-11,12))
#[1.0]
答えは1になりました。実際にこちらで確かめてみましょう。
4次以下の方程式には解の公式があるのでそれを使って解くこともできます。$$ax^3+bx^2+cx+d=0$$について両辺をaで割ると$$x^3+\frac{bx^2}{a}+\frac{cx}{a}+\frac{d}{a}=0$$になります。そして変形すると以下のようになります。$$(x+\frac{b}{3a})^3+(\frac{c}{a}-\frac{b^2}{3a})(x+\frac{b}{3a})+\frac{2b^3}{27a^3}-\frac{bc}{3a^2}+\frac{d}{a}=0$$ここで$$x=y-\frac{b}{3a}$$と置くと$$y^3+(\frac{c}{a}-\frac{b^2}{3a})y+\frac{2b^3}{27a^3}-\frac{bc}{3a^2}+\frac{d}{a}=0$$となります。さらに$$3p=\frac{c}{a}-\frac{b^2}{3a}$$$$2q=\frac{2b^3}{27a^3}-\frac{bc}{3a^2}+\frac{d}{a}$$と置くと$$y^3+3py+2q=0$$となります。これでようやくチルンハウス変換が終わったので次に$$y=u+v$$という一見陳腐な変換を行うと$$y^3+3py+2q=0$$$$(u+v)^3+3p(u+v)+2q=0$$$$u^3+v^3+2q+3(p+uv)(u+v)=0$$となります。そして$$u^3+v^3+2q=0$$$$p+uv=0$$から以下の2式を得ます。$$u^3+v^3=-2q$$$$u^3v^3=-p^3$$するとこれは解と係数の関係から$$t^2+2qt-p^3=0$$の解であることが分かります。また、\(u^3\)と\(v^3\)が対象であることを利用し二次方程式の解の公式より以下の解を得ます。$$u^3=-q+\sqrt{q^2+p^3}$$$$v^3=-q-\sqrt{q^2+p^3}$$これより\(u\)と\(v\)はそれぞれ$$(u_1,u_2,u_3)$$$$=(\sqrt[3]{-q+\sqrt{q^2+p^3}},~\sqrt[3]{-q+\sqrt{q^2+p^3}}w,~\sqrt[3]{-q+\sqrt{q^2+p^3}}w^2)$$$$(v_1,v_2,v_3)$$$$=(\sqrt[3]{-q-\sqrt{q^2+p^3}},~\sqrt[3]{-q-\sqrt{q^2+p^3}}w,~\sqrt[3]{-q-\sqrt{q^2+p^3}}w^2)$$と表せます。これらを\(y=u+v\)に代入し、さらに\(x=y-\frac{b}{3a}\)に代入して得られた解のうち\(uv=-p\)を満たすものになります。
(ちなみに\(w=\frac{-1\pm\sqrt{3}i}{2}\)です。)
なお、ヨビノリさんの動画に三次方程式の解の公式に関するわかりやすい解説動画があったので参考にしてみてください。
それでは実際にコードを書いていきます。まずはp,q,t,wを定義していきます。
def cubic_cardano(a,b,c,d):
p = -b**2/(9.0*a**2) + c/(3.0*a)
q = b**3/(27.0*a**3) - b*c/(6.0*a**2) + d/(2.0*a)
t = complex(q**2+p**3)
w =(-1.0 +1j*3.0**0.5)/2.0ここで最後の行の1jのjは虚数単位を表しています。次にuとvの解を求めていきます。ここではiが0～2のときそれぞれwをi乗するというfor文で書いています。
u = [0,0,0]
v = [0,0,0]
for i in range(3):
u[i] = (-q +t**0.5)**(1.0/3.0) * w**i
v[i] = (-q -t**0.5)**(1.0/3.0) * w**i最後に求めた解から条件に合う解のみを抽出していきます。なお\(uv+p=0\)だと条件が厳しすぎるので0.00001以下になるような条件分岐にしています。そして今回のテーマである実数解を取り出すために虚部が限りなく0に近い場合にそれを実数解とみなしてリストに格納していきます。
x_list = []
for i in range(3):
for j in range(3):
if abs(u[i]*v[j] + p) < 0.00001:
y = u[i] + v[j]
if abs(y.imag) < 0.0000001:
x = y.real - b/(3.0*a)
x_list.append(x)y.imagとy.realはそれぞれyの虚部と実部を表しています。これでx_listというリストに実数解が格納されました。全体のコードは以下のようになっています。
def cubic_cardano(a,b,c,d):
p = -b**2/(9.0*a**2) + c/(3.0*a)
q = b**3/(27.0*a**3) - b*c/(6.0*a**2) + d/(2.0*a)
t = complex(q**2+p**3)
w =(-1.0 +1j*3.0**0.5)/2.0
u = [0,0,0]
v = [0,0,0]
for i in range(3):
u[i] = (-q +t**0.5)**(1.0/3.0) * w**i
v[i] = (-q -t**0.5)**(1.0/3.0) * w**i
x_list = []
for i in range(3):
for j in range(3):
if abs(u[i]*v[j] + p) < 0.00001:
y = u[i] + v[j]
if abs(y.imag) < 0.0000001:
x = y.real - b/(3.0*a)
x_list.append(x)
return x_list
if __name__ == '__main__':
print(cubic_equation(1,-3,7,-5)) #[1.0000000000000002]
DKA法とはデュラン＝カーナー法のことでNewton法に似たようなアルゴリズムでn次方程式の解を求めることができます。このサイトを参考に少しだけカスタマイズしました。
import numpy as np
def cubic_dka(vec,is_complex=False):
dim =len(vec)
if is_complex:
A = np.zeros((dim,dim),dtype=complex)
else:
A = np.zeros((dim,dim))
A[np.arange(dim-1),1+np.arange(dim-1)] =1
A[-1,:] = -vec
ans,vec = np.linalg.eig(A)
return ansここで上記関数は$$x^n+a_{n-1}x^{n-1}+・・・+a_1x+a_0=0$$について各次数の係数を格納しているリストvec[ \(a_0\), \(a_1\), ・・\(a_{n-1}\)]を渡して解いているのでこれを一般的な$$a_nx^n+a_{n-1}x^{n-1}+・・・+a_1x+a_0=0$$についてvec[ \(a_n\), \(a_{n-1}\),・・ \(a_1\), \(a_0\)]を渡して解くように修正したコードが以下になります。
def cubic_dka(vec,is_complex=False):
dim =len(vec)
if is_complex:
A = np.zeros((dim,dim),dtype=complex)
else:
A = np.zeros((dim,dim))
A[np.arange(dim-1),1+np.arange(dim-1)] =1
A[-1,:] = -vec
ans,vec = np.linalg.eig(A)
x_list = []
#実数解を得る
for i in ans:
if np.iscomplex(i) == False or np.iscomplex(i) and abs(i.imag) < 0.0000001:
x_list.append(i.real)
return x_list
#x^3−3x^2+7x−5=0
v = [1,-3,7,-5]
#n次の係数は1にする
div = v.pop(0)
v = [n/div for n in v]
v.reverse()
vec0 =np.array(v)
if __name__ == '__main__':
print(solve(vec0)) #[1.0000000000000002]この手法は何次の方程式でもある程度の精度で解を求めることができるので知っておいて損はない手法だと言えます。
実はn次方程式は1行で解けるんです！Pythonの数値計算ライブラリのnumpyを使って数値的に求めることができます。実際に$$x^3-3x^2+7x-5=0$$について解を求めてみましょう。
import numpy as np
np.roots([1,-3,7,-5])他のどの手法よりも1番短いすっきりとしたコードなのでとりあえず方程式の数値解が知りたいというときはいいかもしれませんね。
また、実数解を求めるにはSympy編と同様の手順で条件分岐します。
x_list = []
for i in x:
if np.iscomplex(i) == False or np.iscomplex(i) and abs(i.imag) < 0.0000001:
x_list.append(i.real)
print(x_list)
#[0.9999999999999997]他の手法と違うところは1未満の数値解が出力されるところでしょうか。おそらく近似法の違いによる誤差なのだと思います。
実際にこれまでご紹介した手法による計算処理の時間を計ってみます。処理時間の計測方法についてはPythonで計算処理時間を計測する方法にまとめましたので詳しく知りたい方は参考にしてみてください。
今回は関数を渡して時間を計測する関数citimeを作りました。また、誤差を無視するために各手法のn回平均の時間を求めています。実際に書いたコードは以下になります。
def ctime(f,n):
t = 0
for _ in range(n):
st = time.perf_counter()
f(1,-3,7,-5)
t += time.perf_counter() - st
res = t / n
return f"{f.__name__} : {res:.10f}"時間計測には無難にtime.perf_counterを用いました。ジュピターノートブックであれば%%timeでもよいのですが余計な処理を省いた純粋な処理時間を計測したかったのでtimeライブラリを用いました。
今回は100回平均の比較を行いました。以下が結果です。
| 手法 | 計測時間 |
| Sympy | 3.21E-03 |
| カルダノ | 7.82E-06 |
| DKA | 5.37E-05 |
| Numpy | 8.90E-05 |
この結果から一番早いのはカルダノの公式を用いた手法でその次にDKA法、その次がnumpyを用いた手法、最後がsympyを用いた手法であることが分かります。もちろん係数により結果が変わることはあるかとは思いますがカルダノの手法が他の手法よりも10乗のオーダーレベルで小さいことから一番処理速度が速い解法であることは揺るがないように思えます。
以下は全体のコードになります。（実数解判定のところは一部関数にまとめました）
import time
from sympy import *
import numpy as np
def cubic_sympy(a,b,c,d):
x = symbols("x")
ans = solve(Eq(a * x ** 3 + b * x ** 2 + c * x + d, 0), x)
x_list = []
for i in ans:
tmp = complex(i)#sympyのAdd型をcomplex型に直してimag分岐を有効化させるため
#実数解だけを抽出
if np.iscomplex(tmp) == False or np.iscomplex(i) and abs(tmp.imag) < 0.0000001:
x_list.append(tmp.real)
return x_list
def cubic_cardano(a,b,c,d):
p = -b**2/(9.0*a**2) + c/(3.0*a)
q = b**3/(27.0*a**3) - b*c/(6.0*a**2) + d/(2.0*a)
t = complex(q**2+p**3)
w =(-1.0 +1j*3.0**0.5)/2.0
u = [0,0,0]
v = [0,0,0]
for i in range(3):
u[i] = (-q +t**0.5)**(1.0/3.0) * w**i
v[i] = (-q -t**0.5)**(1.0/3.0) * w**i
x_list = []
for i in range(3):
for j in range(3):
if abs(u[i]*v[j] + p) < 0.0001:
y = u[i] + v[j]
if abs(y.imag) < 0.0000001:
x = y.real - b/(3.0*a)
x_list.append(x)
return x_list
def cubic_dka(a,b,c,d,is_complex=False):
v = [a,b,c,d]
#n次の係数は1にする
div = v.pop(0)
v = [n/div for n in v]
v.reverse()
vec =np.array(v)
dim =len(vec)
if is_complex:
A = np.zeros((dim,dim),dtype=complex)
else:
A = np.zeros((dim,dim))
A[np.arange(dim-1),1+np.arange(dim-1)] =1
A[-1,:] = -vec
ans,vec = np.linalg.eig(A)
#実数解を得る
return check_rns(ans)
def cubic_numpy(a,b,c,d):
x = np.roots([a,b,c,d])
return check_rns(x)
def check_rns(x):
x_list = []
for i in x:
if np.iscomplex(i) == False or np.iscomplex(i) and abs(i.imag) < 0.0000001:
x_list.append(i.real)
return x_list
def ctime(f,n):
t = 0
for _ in range(n):
st = time.perf_counter()
f(1,-3,7,-5)
t += time.perf_counter() - st
res = t / n
return f"{f.__name__} : {res:.10f}"
if __name__ == '__main__':
n = 100
print(ctime(cubic_sympy,n))
print(ctime(cubic_cardano,n))
print(ctime(cubic_dka,n))
print(ctime(cubic_numpy,n))
三次方程式を解く方法は様々ありまが今回の結果からカルダノの公式の手法が最も早く、Sympyを用いる手法が最も遅いということが分かりました！私は研究で三次方程式を大量に解かなければならないのですが、Sympyでは20分かかっていたものを数秒に抑えることに成功しました。なのでもし早い処理が必要な場合はカルダノの公式による手法をお勧めします！
もしもっと早い手法があるよ！っていう場合はぜひコメントしていただけると嬉しいです。