【Python】多項式代数方程式の解を一気に見つけるアルゴリズムの実装【数値計算】
はじめに
学校のレポートをやるために、$ z^{ n } + a_{ n - 1 } z^{ n - 1 } + \cdots + a_0 = 0 \, ( z \in \mathbb{ C } ) $という多項式代数方程式の全ての根を数値的に求める必要が出たのでPythonで実装しました。
Newton法とか、1個だけ求めるアルゴリズムは知っていはいましたが、全て一気に解ける方法はないものかと調べていたら...、あったあった。先日買った「数値計算の常識」という本に解き方が載っていました。この本いいですね。気になる人の為にAmazonのリンク貼っておきます。(「数値計算の常識」のAmazonリンク)
なお、Pythonの数値計算・数値シミュレーションの優秀なライブラリであるScipyにも代数方程式をとくメソッドがありましたが、今回はNumpyのみを使って自作してみました。
調べてたら、Numpyにもありました。多項式代数方程式を解く方法。Python優秀すぎる。まじでライブラリ使わないで実装する意味のなさ...。
という訳で、Scipy、Numpyのライブラリを用いて実装する方法については、おまけに記しました。
解きたい方程式
解きたい方程式を確認しておきます。次の式(1)です
$$
p(z) = z^{ n } + a_{ n - 1 } z^{ n - 1 } + \cdots + a_0 = 0 \hspace{0.3cm} ( z \in \mathbb{ C } ) \tag{1}
$$
連立方程式を解くアルゴリズム
今から記すアルゴリズムは連立法と言われるものらしいです。
では、連立法の流れをみていきましょう。
連立法の流れ
連立法は、基本、反復法のアルゴリズムまんまです。
すなわち、$n$個の$\nu$次の近似解$z_1^{(\nu)}, ..., z_n^{(\nu)}$に対して、ある$n$次関数$\varphi_1(z_1,...,z_n),...,\varphi_n(z_1,...,z_n)$を用いて$\nu+1$次の近似解
$$
z_i^{(\nu +1)} = z_i^{(\nu)} + \varphi_i(z_1^{(\nu)}, ..., z_n^{(\nu)}) \tag{2}
$$
で計算して、ある正数$\varepsilon$について、ある$i$が$|p(z_i^{(\nu})| < \varepsilon$となったら反復を停止するものです。
これには、以下の手順が含まれます。
- 初期値$(z_1^{(0)}, ..., z_n^{(0)})$を定める
- 反復関数$\varphi_i(z_1, ..., z_n)$の定めて、$\nu$近似解$(z_1^{(\nu)}, ..., z_n^{(\nu)})$から$\nu + 1$近似解$(z_1^{(\nu+1)}, ..., z_n^{(\nu+1)})$を作る。
- 作られた$\nu + 1$近似解に対して誤差を見積もり反復を停止するか否かを判定する
初期値$(z_1^{(0)}, ..., z_n^{(0)})$を定める
では、まずどう初期値を求めるのがいいでしょうか。
初期値は求める値、すなわち、方程式の根から近い方がいいです。そのため、根が存在する範囲をまず調べます。
そのために、根の重心$b$を求め、その重心から全ての根を含む最大半径$r$を得る必要があります。
根の重心は、解と係数の関係から、解を $ \alpha_{ 1 }, \alpha_{ 2 }, ..., \alpha_{ n }$ とすると、
$$
b = - \frac{ a_{ n - 1 } }{ n } = - \frac{ \alpha_{ 1 } + \alpha_{ 1 } + \cdots + \alpha_{ n } }{ n } \tag{3}
$$
で求められます。
半径$r$の求めかたはテクニカルです。
すなわち、根の重心周りで$p(z) $を展開した式
$$
p(\zeta + b) = \zeta^n + c_{n-1} \zeta^{n-1} + \cdots + c_1 \zeta + c_0 \tag{4.1}
$$
の係数を組立除法(あとで説明します。)(もしくはテイラー展開？)で求めて、実係数の方程式
$$
x^n - |c_{n-2}| x^{n-2} - \cdots - |c_1|x - |c_0| = 0 \tag{4.2}
$$
を考えると、式(4)を$x^n$でくくると、
$$
x^n \{ 1 - ( \frac{|c_{n-2}|}{x^2} +\cdots + \frac{|c_1|}{x^{n-1}} + \frac{|c_0|}{x^n} ) \} = 0
$$
となり、
$$
f(x) = \frac{|c_{n-2}|}{x^2} +\cdots + \frac{|c_1|}{x^{n-1}} + \frac{|c_0|}{x^n}
$$
は$x>0$で$f(x) > 0$でかつ単調減少であり、$x<0$で$f(x) < 0$でかつ単調増加ゆえ、$g(x) = 1$と$x = r > 0$で実数の範囲内でただ一つの共有点を持つことがわかります。ゆえに、式(4)にはただ一つの正の解があり、それを$r$とします。
そして、根の存在範囲を得られたら、次の初期値から出発するのがいいらしいです。(なんでかは知りません。調べてないです。すいません。)
$$
z_i^{0} := b + r \exp [j (\frac{2\pi (i-1)}{n} + \frac{3}{2n})] \tag{5}
$$
ただし、$j$は虚数単位です。
反復関数の定め方
次に反復関数$\varphi_i(z_1, ..., z_n)$をいかに定めるかです。これこそ連立法のキモです。
今回は、連立法でも2次法と言われるものを用います。これは、Newton法の反復法の2次収束の考え方を背景に持ちます。
すなわち、Newton法での反復関数 $$ \varphi_i(z_1, ..., z_n) = - \frac{p(z_i^{(\nu)})}{p'(z_i^{(\nu)})} $$ において、$p(z) = \prod_{k =1 }^n (z- \alpha_k )$とすると、 $$ p'(z) = \sum_{k=1}^n \prod_{j \neq k} (z- \alpha_j ) $$ となるため、$z_i^{(\nu)}$が$\alpha_i$に十分近いとして、$z_i^{(\nu)}- \alpha_i \simeq 0$を用いると、 $$ p'(z) = \prod_{j \neq i} (z_i^{(\nu)}- \alpha_i) $$ となり、他の$j \neq i$についても$z_j^{(\nu)} \simeq \alpha_j$を用いると、 $$ p'(z) = \prod_{j \neq i} (z_i^{(\nu)}- z_j^{(\nu)}) $$ となるため、結局、反復関数を以下のようにおきます。 $$ \varphi_i(z_1, ..., z_n) := - \frac{p(z_i)}{\prod_{k = 1, \, k \neq i}^n (z_i - z_k)} \tag{6} $$
今回は、連立法でも2次法と言われるものを用います。これは、Newton法の反復法の2次収束の考え方を背景に持ちます。
すなわち、Newton法での反復関数 $$ \varphi_i(z_1, ..., z_n) = - \frac{p(z_i^{(\nu)})}{p'(z_i^{(\nu)})} $$ において、$p(z) = \prod_{k =1 }^n (z- \alpha_k )$とすると、 $$ p'(z) = \sum_{k=1}^n \prod_{j \neq k} (z- \alpha_j ) $$ となるため、$z_i^{(\nu)}$が$\alpha_i$に十分近いとして、$z_i^{(\nu)}- \alpha_i \simeq 0$を用いると、 $$ p'(z) = \prod_{j \neq i} (z_i^{(\nu)}- \alpha_i) $$ となり、他の$j \neq i$についても$z_j^{(\nu)} \simeq \alpha_j$を用いると、 $$ p'(z) = \prod_{j \neq i} (z_i^{(\nu)}- z_j^{(\nu)}) $$ となるため、結局、反復関数を以下のようにおきます。 $$ \varphi_i(z_1, ..., z_n) := - \frac{p(z_i)}{\prod_{k = 1, \, k \neq i}^n (z_i - z_k)} \tag{6} $$
こうすると、反復列は2次収束します。
停止条件 : 誤差を見積もる(Smithの方法)
そうして作った反復列に対し、停止条件を考えるために、誤差の見積もりをします。
式(1)の解の近似値$z_1,...,z_n$があるとき、真の存在範囲を厳密に求める方法があります。この方法の根拠は長くなるので割愛しますが、大まかに次の方法で求められます。
まず、全ての$i$について、$|p(z_i)| < \delta_i$となるような$\delta_i$を計算して、それを用いて$\varphi_i$の上界$\tilde{\varphi_i}$を $$ \varphi_i(z_1,...,z_n) \le \tilde{\varphi_i}(z_1,...,z_n) := \frac{|p(z_i)| + \delta_i}{\prod_{j \neq i} |z_i - z_j|} $$ と算出しておきます。
そうして、$n$個の円、 $$ C_i := \{ z | |z_i + \varphi_i (z_1,...,z_n) - z| \le (n-1) \tilde{\varphi_i}(z_1,...,z_n) \} $$ を作ると、全ての根は$\cup_{i=1}^n C_i$の中にあります。
この$C_i$の半径$(n-1) \tilde{\varphi_i}(z_1,...,z_n)$が十分に小さい正数$\varepsilon$を用いて、 $$ (n-1) \tilde{\varphi_i}(z_1,...,z_n) < \varepsilon $$ となった時に、反復列は収束したとみなします。
まず、全ての$i$について、$|p(z_i)| < \delta_i$となるような$\delta_i$を計算して、それを用いて$\varphi_i$の上界$\tilde{\varphi_i}$を $$ \varphi_i(z_1,...,z_n) \le \tilde{\varphi_i}(z_1,...,z_n) := \frac{|p(z_i)| + \delta_i}{\prod_{j \neq i} |z_i - z_j|} $$ と算出しておきます。
そうして、$n$個の円、 $$ C_i := \{ z | |z_i + \varphi_i (z_1,...,z_n) - z| \le (n-1) \tilde{\varphi_i}(z_1,...,z_n) \} $$ を作ると、全ての根は$\cup_{i=1}^n C_i$の中にあります。
この$C_i$の半径$(n-1) \tilde{\varphi_i}(z_1,...,z_n)$が十分に小さい正数$\varepsilon$を用いて、 $$ (n-1) \tilde{\varphi_i}(z_1,...,z_n) < \varepsilon $$ となった時に、反復列は収束したとみなします。
ごちゃごちゃ言いましたが、実際には、$|p(z_i^{(\mu)})|$の計算をして、それが、十分に小さい正数$\varepsilon$を用いて、
$$
\forall i \in \{1,2,..., n\} \hspace{0.3cm} |p(z_i^{(\mu)})| < \varepsilon
$$
となる時に反復を停止する方法でもいいです。
以上が連立法の流れとなります。
(補足) 組み立て除法
ここで式(2)を得るために用いる組立除法について軽く触れておきます。
組み立て除法とは、テーラー展開した時の係数を求めるアルゴリズムといっても過言ではなさそうです。
すなわち、
$$
\begin{array}{ccl}
p(x) &=& p(\alpha) + p'(\alpha)(x-\alpha) + \frac{p''(\alpha)}{2!}(x-\alpha)^2 + \cdots tag{7} \\\
&=& c_0 + c_1 \zeta + c_2 \zeta^2 + \cdots
\end{array}
$$
に対して、係数$c_k = \frac{p^{k}(\alpha)}{k!}$を求めるのに、わざわざ微分してべき乗計算するのは計算コストが高いので、工夫をしましょう。式(7)より
$$
\begin{array}{ccl}
p(x) &=& p(\alpha) + (x-\alpha) ( p'(\alpha) +\frac{p''(\alpha)}{2!}(x-\alpha) + \cdots ) \\\
&=& p(\alpha) + q^{(1)}(x, \alpha)
\end{array}
$$
となるため、$p(\alpha)$は$p(x)$を$x-\alpha$で割った時の余りとなることがわかります。そのため、$p(x)$を$x-\alpha$で割って余りを求めることを考えましょう。この時、次の表を用いて求められます。(証明は割り算の筆算を書いてみるとわかります)
これにより、$p(\alpha)$は求められました。
また、$q^{(1)}(x, \alpha) = a_0 + b_1^{(1)}x^{n-2} + \cdots + b_{n-1}^{(1)}$となるため、この$q^{(1)}(x, \alpha)$に対して、再び、 $$ q^{(1)}(x, \alpha) = p'(\alpha) + (x-\alpha) q^{(2)}(x, \alpha) $$ とすると、$p'(\alpha)$は$q^{(1)}(x, \alpha)$を$(x-\alpha)$で割った時の剰余で、$q^{(2)}(x, \alpha)$は商。先ほどと同じようにして、それぞれ求められる。これを繰り返して行くことで、 $$ \frac{p^{k}(\alpha)}{k!} = b_{n-k}^{(k+1)} $$ として求められます。
| 係数 | $a_0 $ | $a_1$ | $a_3$ | $\ldots$ | $a_{n-1}$ | $a_{n}$ |
|---|---|---|---|---|---|---|
| $\times \alpha$ | $a_0 \alpha$ | $b_1^{(1)} \alpha$ | $\ldots$ | $b_{n-2}^{(1)} \alpha$ | $b_{n-1}^{(1)} \alpha$ | |
| 和 | $a_0 $ | $b_1^{(1)}$ | $b_2^{(1)}$ | $\ldots$ | $b_{n-1}^{(1)} $ | $b_{n}^{(1)} = p(\alpha)$ |
また、$q^{(1)}(x, \alpha) = a_0 + b_1^{(1)}x^{n-2} + \cdots + b_{n-1}^{(1)}$となるため、この$q^{(1)}(x, \alpha)$に対して、再び、 $$ q^{(1)}(x, \alpha) = p'(\alpha) + (x-\alpha) q^{(2)}(x, \alpha) $$ とすると、$p'(\alpha)$は$q^{(1)}(x, \alpha)$を$(x-\alpha)$で割った時の剰余で、$q^{(2)}(x, \alpha)$は商。先ほどと同じようにして、それぞれ求められる。これを繰り返して行くことで、 $$ \frac{p^{k}(\alpha)}{k!} = b_{n-k}^{(k+1)} $$ として求められます。
実装
という訳で、2次の連立法で多項式方程式を求める場合、実装を次のような流れで進めます。
- 係数配列$a = [a_0, a_1,...,a_{n-1}]$で多項式を定義する
- Newton法の実装する
根の存在範囲を求める。
- 根の重心を求める
- 組み立て除法により、式(4.1)の係数$c_i$を求める
- 式(4.2)の解をNewton法により求める
初期値を求めて反復法を開始して、根の配列$z = [z_0, z_1,..., z_{ n -1 }]$を求める。
このように進めて行きます。
配列coefficientsを係数にもつ多項式関数を作る
def getGfunc(coefficients):
"""N次式を作る
coefficients = [a_0, a_1,...., a_{N-1}]として、a_0 + a_1 x^1 + ... + a_{N-1} x^{N-1}を得る。
--Input---
coefficients : ndarray. N個の要素をもつ
-- return --
g : coefficientsを係数とするzに関するN次式。
"""
def g(x):
func_val = 0
for index in range(len(coefficients)):
func_val += coefficients[index] * (x**(index))
return func_val
return g
複素数版Newton法の実装
Newton法で解を見つけるときは以下のコードでできます。複素数版Newton法です。
def newtonMethod(func, init_val, max_step = 50):
delta_val = complex(0.0001, 0.0001) # 微小区間の幅
# 以下Newton法
for i in range(max_step): # 回数制限をかける
f = func(init_val)
df = (func(init_val + delta_val) - func(init_val - delta_val)) / (2*delta_val) # 微分値の取得
init_val = init_val - f/df # 解の更新
ans_x = init_val # 近似解
ans_y = func(ans_x) # そのときのyの値(確認用)
print('x : '+str(ans_x))
print('y : '+str(ans_y))
return ans_x
def func_y(x_dat):
func = 2*(x_dat**2) - 3*x_dat - 4
return func
if __name__ == '__main__':
init_val = complex(4, 0) # 初期値
print(newtonMethod(func_y, init_val))
最大半径を求める
組立除法で、式(2)の係数$c_i$を求め、さらに、式(2)の正の根$r$をNewton法で求めることを目指す。
def getRadiusBykumitateZyoho(coefficients):
"""組立除法とニュートン法を用いて最大半径を求める
--Input---
coefficients : ndarray. N個の要素をもつ係数配列.
-- return --
g : float. coefficientsを係数とするzに関するN次式の解に対して、
解の重心からの全ての解が含まれる最大半径rを求める。
"""
dim_n = len(coefficients) - 1
# 根の重心
root_centroid = coefficients[-2] / dim_n
# rを求める際にできる関数の係数を格納するリスト
new_coefficients = []
b = 0
for k in range(dim_n):
for i in range(dim_n - k):
b = coefficients[i] + b * root_centroid
new_coefficients.append(-abs(b))
new_coefficients.append(1)
# 新しい関数を作る
new_func = getGfunc(new_coefficients)
radius = newtonMethod(new_func, 10)
return abs(radius)
連立法
# zについてのN次の代数方程式をとく
def solveEquationNdim(coefficients, max_step=10**2):
"""連立法を用いて、代数方程式をとく"""
dim_n = len(coefficients) - 1
root_centroid = - coefficients[-2]/dim_n
# 求める関数
objectFunc = getGfunc(coefficients)
# 最大半径を求める
radius = getRadiusBykumitateZyoho(coefficients)
# 初期値を設定する
z_list = []
for i in range(dim_n):
z_list.append(root_centroid + radius*np.exp(1j *(2*np.pi*((i+1) - 1)/dim_n + 3/(2*dim_n))))
def timesFunc(val_list, index):
# 分母
denominator = 1
val_i = val_list[index]
for j in range(len(val_list)):
if j == index:
continue
denominator *= (val_i - val_list[j])
return denominator
# 求まった根の個数
ans_num = 0
# 更新された根を格納するリスト
new_z_list = [0]*dim_n
# 関数値を格納するリスト
func_val_list = [0]*dim_n
for i in range(max_step):
for i in range(len(z_list)):
fun_val = objectFunc(z_list[i])
func_val_list[i] = fun_val
denominator = timesFunc(z_list, i) # 分母を計算
psi = - fun_val/denominator # 更新する値を計算
new_z = z_list[i] + psi # 更新
new_z_list[i] = new_z
new_radius = new_z + abs(fun_val)/abs(denominator)
z_list = new_z_list
if (np.max(abs(np.array(func_val_list))) < 1e-8):
print(i+1,"回で停止")
return np.sort(np.array(z_list))
print("収束しなかった")
return np.sort(np.array(z_list))
おまけ
Scipy、Numpyのライブラリを用いて実装する方法について軽くみていきます。
Numpyのライブラリで解く方法
Numpyには多項式を扱うためのライブラリがあります。それを使うことで一瞬で解が求まるみたいです。すげえ。
Numpyには多項式を扱う際に、二つの方法があるみたいです。
1個目の方法 : poly1dメソッド
1個は、poly1dメソッドを用いる方法。以下のコードでできます。
以下のコードでは$x^{ 3 } - 2 x^{ 2 } + x - 2 = 0$を解いてます。
# Define a polynomial x^3 - 2 x^2 + x - 2.
p = np.poly1d([1, -2, 1, -2])
# print polynomial
print(p)
# Solve f(x) == 0.
print("roots: ", p.roots)
こうすると、
3 2 1 x - 2 x + 1 x - 2 roots: [ 2.00000000e+00+0.j -4.16333634e-16+1.j -4.16333634e-16-1.j]
と出力されます。
ちなみに、p(0)とかすると、$x=0$を代入した値を得られるみたい。
2個目の方法 : polynomialメソッドを使う
polynomialメソッドもべき乗多項式を扱うものです。実装例は以下のよう。
なお、先ほどと同じく、$x^{ 3 } - 2 x^{ 2 } + x - 2 = 0$を解いてます。
from numpy.polynomial.polynomial import Polynomial
# Define a polynomial x^3 - 2 x^2 + x - 2.
testEqF = Polynomial([-2, 1, -2, 1])
print("polynomial: ", testEqF)
# Solve f(x) == 0.
print("roots: ", testEqF.roots())
注意すべきは、先ほどと係数の指定する順番が真逆だということ。
出力は、
polynomial: poly([-2. 1. -2. 1.]) roots: [-1.83880688e-16-1.j -1.83880688e-16+1.j 2.00000000e+00+0.j]
となりました。根の表示は先ほどとは少し順番が異なりますね。
Scipyのライブラリでとく方法
こちらはScipyのシンボリック機能とsolveメソッドを使う方法。
お次は、$2x^{ 3 } + 3x - 5 = 0$を解いてみます。以下のコードででできます。
from sympy import *
x=Symbol('x') # 文字'x'を変数xとして定義
"""
方程式: solve
例: 2 * x **3 +3*x-5=0の根を求める
"""
sol=solve(2 * x **3 +3*x-5, x)
print(sol)
出力は、
[1, -1/2 - 3*I/2, -1/2 + 3*I/2]
Iというのは虚数を表すのでしょう。こっちの方が直感的。
最後に
今回の記事、思ったより長かった...。疲れた...。
まあ、でもいつかまた必要になると思うからね。。。。
でも、まじでライブラリ使わないならPython使う意味ほぼないし、C++あたりで実装し直そうかな...。
全く関係ないけど、この記事を書きながらスマブラやってました。楽しい。
キングクルール使いやすくて良いね。
ということで、今回はここまでです。バイバーイ。
参考文献
連立法
Numpyライブラリで解く
Scipyで代数方程式を解く
