【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)です
連立方程式を解くアルゴリズム
今から記すアルゴリズムは連立法と言われるものらしいです。
では、連立法の流れをみていきましょう。
連立法の流れ
連立法は、基本、反復法のアルゴリズムまんまです。
これには、以下の手順が含まれます。
- 初期値$(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$を得る必要があります。
半径$r$の求めかたはテクニカルです。
そして、根の存在範囲を得られたら、次の初期値から出発するのがいいらしいです。(なんでかは知りません。調べてないです。すいません。)
反復関数の定め方
今回は、連立法でも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の方法)
そうして作った反復列に対し、停止条件を考えるために、誤差の見積もりをします。
まず、全ての$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 $$ となった時に、反復列は収束したとみなします。
以上が連立法の流れとなります。
(補足) 組み立て除法
ここで式(2)を得るために用いる組立除法について軽く触れておきます。
組み立て除法とは、テーラー展開した時の係数を求めるアルゴリズムといっても過言ではなさそうです。
係数 | $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))
最大半径を求める
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で代数方程式を解く