n次方程式って
x^n+a_{n-1}x^{n-1}+\cdots+a_1x+a_0 =0
を解きたいシチュエーションに出会った。これは一般にn次方程式と呼ばれている気がする。（1変数の代数方程式ともいうのか）
これは5次以上のものについては代数的な解がないことが知られている。
とはいっても5次についてはいろんな特殊関数を使えば表現できるようだし、近似解については数値計算もできる。
4次以上についてはそもそも数値計算を用いるのが普通のようだ。
Matlabなんかではrootsという代数方程式の数値計算を行う関数が用意されているようだが、pythonでは特にライブラリが見当たらなかった。というわけで実装してみることにした。
アルゴリズム
解を求めるとなるとNewton法が最初に思いつくように思われるが、n次方程式の解を求めるNewton法に近いアルゴリズムとしてDKA法というものが知られているようだ。ただこのアルゴリズムを実装するのは面倒だし、Matlabの実装をみると、固有値計算に帰着する方法を用いているらしい。今回はこれを用いることにした。
固有値を定義に沿って計算するときは
\mathrm{det}(sI-A)=0
を解くことになるわけだが、これは結局行列がn次の時n次方程式となる。固有値計算については様々な方法があるようで、numpyにおいてもnp.linalg.eigなる関数が用意されている。よって解を求めたい方程式を固有多項式としてもつ行列が分かればいいわけだが、これについてはコンパニオン行列という以下のようなものが知られている。
A=
\left(
\begin{array}{ccccc}
0 & 1 & 0&\ldots&0 \\
\vdots & \ddots & 1&\ddots&\vdots \\
\vdots & \ldots & \ddots&\ddots&0\\
0&\ldots&\ldots&0&1\\
-a_0&-a_1&\ldots&-a_{n-2}&-a_{n-1}
\end{array}
\right)
これの固有方程式は
|sI-A| = s^n+a_{n-1}s^{n-1}+a_{n-2}s^{n-2}+\ldots+a_1s+a_0
となる。
よって任意係数をもつn次方程式について以上のコンパニオン行列の固有値をすべて計算してやれば、解が求まることになる。
なお最初に紹介したDKA法についてはこちらを参照。実装も落ちていた。
実装
import numpy as np
def solve(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
vec0 =np.array([-120,274,-225,85,-15])
vec1 =np.array([1,0])
vec2 =np.array([-1,5,-10,10,-5])
print(solve(vec0))
print(solve(vec1))
print(solve(vec2))
[ 5. 4. 3. 2. 1.]
[ 0.+1.j 0.-1.j]
[ 1.00079742+0.00057977j 1.00079742-0.00057977j 0.99969502+0.00093688j
0.99969502-0.00093688j 0.99901512+0.j ]
係数が虚数を含む場合は、is_complexをTrueにしないと、勝手に係数が実数に変換されてしまうことに注意。なお常にdtype=complexとしていると、実数固有値が実数固有値として計算されにくくなるようなので、普段はfloatで計算している(dtypeによってアルゴリズムが違うのだろうか？）
例はそれぞれ
(x-1)(x-2)(x-3)(x-4)(x-5)=x^5-15x^4+85x^3-225x^2+274x-120=0\\
x^2+1=(x-i)(x+i)=0\\
(x-1)^5=0
を計算している。解の公式のない5次方程式についても重解含めそれなりの精度で計算されていることが分かる。