Jenkins-Traubアルゴリズムで多項式の根を求める～Pythonコード付き～　　

多項式の根を求めるアルゴリズムの一つ、Jenkins-Traub法(ジェンキンズ-トラウブ法)の導出方法と、Pythonでの実行コードについて解説する。

1. Jenkins-Traubアルゴリズムの導出
2. Jenkins-Traub法のPythonコード

Jenkins-Traubアルゴリズムの導出

求根アルゴリズムの一つであるJenkins-Traub法とは、$n$次の多項式$p(z)$の根、すなわち$p(z)=0$の解$z$を求める手法のひとつである。

$$ p(z) = z^n + a_1 z^{n-1} + \cdots + a_{n-1}z + a_n $$

$p(z)$の根を$\alpha_j$、重根度を$m_j$とすると、$p(z)$は以下のように表せる。

$$p(z) =\prod_{j=1}^k (z – \alpha_j)^{m_j}$$

$p(z)$を$z$で微分する。

$$ \begin{eqnarray}
p'(z) &=& \sum_{j=1}^k m_j \frac{p(z)}{z-\alpha_j} \\
&=& \sum_{j=1}^k m_j P_j(z) \\
\Bigl( P_j(z) &=& \frac{p(z)}{z-\alpha_j} \Bigr) \tag{1} \\
\end{eqnarray} $$

ここで次のような多項式列$ \{H^{(\lambda)}(z) \}_{\lambda} $を考える。

$$ \begin{eqnarray}
H^{(0)}(z) &=& p'(z) \\
H^{(\lambda)}(z) &=& \sum_{j=1}^k c_j^{(\lambda)} P_j (z) \tag{2} \\
\end{eqnarray} $$

係数$c_j^{(\lambda)} $については、下記を満たすようなものをとる。

$$ \begin{eqnarray}
c_j^{(0)} &=& m_j \ \ (j=1,2,\cdots ,k) \\
\Bigl( d_j^{(\lambda)} &=& \Bigr) \frac{c_j^{(\lambda)}}{c_1^{(\lambda)}} \xrightarrow[\lambda \xrightarrow{} \infty]{} 0 \ \ (j=2,\cdots ,k) \tag{3}\\
\end{eqnarray} $$

すると$H^{(\lambda)}(z)$は、

$$ H^{(\lambda)}(z) \xrightarrow[\lambda \xrightarrow{} \infty]{} c_1^{(\lambda)}P_1(z) \tag{4} $$

となる。

このとき、次式で定義される数列$ \{ t_{\lambda} \}_{\lambda} $は根$\alpha_1$に近づく。(＊)

$$ \begin{eqnarray}
t_{\lambda +1} &=& s_{\lambda} \ – \ \frac{p(s_\lambda)}{\widetilde{H}^{(\lambda +1)}(s_{\lambda})} \tag{5} \\
\widetilde{H}^{(\lambda)}(z) &=& \frac{H^{(\lambda)}(z)}{\sum_{j=1}^k c_j^{(\lambda)}} \tag{6}
\end{eqnarray} $$

式(5)の$ \{ s_{\lambda} \}_{\lambda}$は任意の複素数列である。

(＊)の証明

式(5)に式(6)を代入すると、

$$ t_{\lambda +1} = s_{\lambda} \ – \ \frac{ \sum_{j=1}^k c_j^{(\lambda +1)} p(s_\lambda)}{H^{(\lambda +1)}(s_{\lambda})} $$

となる。さらに、式(1)で$j=1$とした$ \ p(s_{\lambda})=(s_{\lambda} – \alpha_1)P_1(s_{\lambda}) \ $と式(2)を代入すると、

$$ \begin{eqnarray}
t_{\lambda +1} &=& s_{\lambda} \ – \ \frac{ \sum_{j=1}^k c_j^{(\lambda +1)} (s_{\lambda} – \alpha_1)P_1(s_{\lambda}) }{ \sum_{j=1}^k c_j^{(\lambda +1)} P_j(s_{\lambda})} \\
&=& s_{\lambda}\ – \ \frac{ (s_{\lambda}-\alpha_1) P_1(s_{\lambda}) c_1^{(\lambda+1)} (1+\sum_{j=2}^k d_j^{(\lambda +1)} ) }{ c_1^{(\lambda +1)} P_1(s_{\lambda}) ( 1+ \sum_{j=2}^k d_j^{(\lambda +1)} \frac{P_j(s_{\lambda})}{P_1(s_{\lambda})} ) } \\
& \xrightarrow[\lambda \xrightarrow{} \infty]{} & s_{\lambda} \ – \ (s_{\lambda} – \alpha_1) = \alpha_1
\end{eqnarray} $$

((＊)の証明終了)

次に$ \widetilde{H}^{(\lambda)} (z)$および$H^{(\lambda)}(z) $の具体的な形について考えてみる。

式(1)より$ P_j(z) $の最大次数は$ (n-1) $であり、その$ (n-1) $次の項の係数は$ 1$である。そして、式(2)より$H^{(\lambda)}(z)$も$ (n-1) $次の多項式である。$ \widetilde{H}^{(\lambda)} (z)$については式(6)より$H^{(\lambda)}(z)$と同様に$ (n-1) $次の多項式であり、さらに式(6)に式(3),(4)を適用することで、

$$ \widetilde{H}^{(\lambda)} (z) \xrightarrow[\lambda \xrightarrow{} \infty]{} P_1(z) $$

となるので、$\lambda$が十分大きいときには$ \widetilde{H}^{(\lambda)} (z)$の$ (n-1) $次の項の係数は$ P_1(z) $と同様に$1$となる。すなわち$ \widetilde{H}^{(\lambda)} (z)$は、$H^{(\lambda)}(z)$を$ (n-1) $次の項の係数が$1$になるように規格化することで求められる。

そして、式(4)を満たすような多項式列$ \{H^{(\lambda)}(z) \}_{\lambda} $については次のような漸化式を用いてやればよい。

$$ H^{(\lambda+1)}(z) = \frac{1}{z-s_{\lambda}} \Bigl[ H^{(\lambda)}(z) – \frac{H^{(\lambda)}(s_{\lambda})}{p(s_{\lambda})}p(z) \Bigr] \tag{7} $$

[・・・]内は$z=s_{\lambda}$の時$0$になるため$ (z-s_{\lambda}) $で割り切れる。すなわち、すべての$ \lambda$について式(7)で定義される多項式列$ \{H^{(\lambda)}(z) \}_{\lambda} $は$ (n-1) $次の多項式である。

式(7)より

$$ \begin{eqnarray}
H^{(\lambda+1)}(z) &=& \frac{p(z)}{z-s_{\lambda}} \Bigl[ \ \sum_{j=1}^{k} \frac{c_j^{(\lambda)}}{z-\alpha_j} \ – \ \sum_{j=1}^{k} \frac{c_j^{(\lambda)} }{s_{\lambda}-\alpha_j} \ \Bigr] \\
&=& \frac{p(z)}{z-s_{\lambda}} \sum_{j=1}^{k} \frac{ c_j^{(\lambda)} (s_{\lambda}-z) }{(z-\alpha_j)(s_{\lambda}-\alpha_j) } \\
&=& p(z)\sum_{j=1}^{k} \frac{ c_j^{(\lambda)} }{(z-\alpha_j)(\alpha_j-s_{\lambda} )} \\
&=& \sum_{j=1}^{k} P_j(z)\frac{ c_j^{(\lambda)} }{(\alpha_j-s_{\lambda} ) } \\
\end{eqnarray} $$

式(2)で$\lambda \xrightarrow{} \lambda+1 $ としたものと上式を比較すると、

$$ \begin{eqnarray}
c_j^{(\lambda+1)} &=& \frac{ c_j^{(\lambda)} }{ \alpha_j \ – \ s_{\lambda} } \\
&=& \frac{ c_j^{(\lambda-1)} }{ (\alpha_j \ – \ s_{\lambda})(\alpha_j \ – \ s_{\lambda-1}) } \\
& \cdots& \\
&=& \frac{c_j^{(0)}}{\prod_{l=0}^{\lambda} (\alpha_j \ – \ s_l)} \\
&=& \frac{m_j}{\prod_{l=0}^{\lambda} (\alpha_j \ – \ s_l )}
\end{eqnarray} $$

$\{s_{\lambda}\}$は任意の複素数列なので、根$\alpha_1$に収束する数列$\{t_{\lambda}\}$とすれば

$$ \begin{eqnarray}
|c_1^{(\lambda)}| &\xrightarrow{}& \infty \ \ (\lambda \xrightarrow{} \infty)\\
|c_j^{(\lambda)}| &<& \infty \ \ (j=2,3,\cdots,k) \\
\end{eqnarray} $$

となり、式(3)および式(4)が満たされていることがわかる。

参考文献：

Mekwi, Wankere R. “Iterative methods for roots of polynomials.” (2001).

Jenkins–Traub algorithm – Wikipedia

関連記事：

ニュートン法による複素求根アルゴリズム~Pythonコード付き~

ラゲール法による求根アルゴリズム～Pythonコード付き～

Jenkins-Traub法のPythonコード

以下のプログラムでは、

$$p(z) = z^5 +1 $$

を例に、初期値$z_0$を$\{ z \mid -10<Re(z)<10, -10<Im(z)<10 \} $の範囲でJenkins-Traub法を適用し、根の偏角の大きさをカラーマップで表示しています。

# -*- coding: utf-8 -*-
"""
『Jenkins Traub法』で複素平面の各点を初期値としてp(z)の根を求め、その偏角をカラーマップにする
"""
import numpy as np
import matplotlib.pyplot as plt

LIMIT = 100  # 計算回数の上限 
EPS = 1E-10  # ゼロ点の判定範囲. abs(p(z)) < eps ならばゼロ点とする.

p = np.poly1d([1,0,0,0,0,1])  # p(z) = z^5 + 1
dp = p.deriv()  # 微分
print(p)

# 初期値の計算範囲
real_min=-10
real_max=-real_min
im_min=-10
im_max=-im_min
# ステップ数
N_real=500
N_im=500


# 複素平面の各点を初期値として、p(z)の根を求めその偏角をカラーマップにする
def plot_complex(complex_func):
    x, y = np.meshgrid(np.linspace(real_min,real_max, N_real), np.linspace(im_min,im_max, N_im))
    z = x + y*1j
    
    solutions = np.zeros_like(z)
    for ix in range(N_real):
        for iy in range(N_im):
            solutions[iy][ix] = Root_of_func(z[iy][ix])
        print(ix,"/",N_real-1)
            
    angles = (np.angle(solutions) + 2*np.pi) % (2*np.pi) *180/np.pi 
    plt.figure(figsize=(4, 4), dpi=300)
    im = plt.imshow(angles, origin="lower",cmap="hsv",vmin=0,vmax=360,\
                    extent=[real_min,real_max,im_min,im_max]) 
    plt.colorbar(im, fraction=0.046, pad=0.04) 
    return angles


# Jenkins Traub法
def Root_of_func(z):
    H = dp  # 多項式H0(z)
    s = z  # s0
    
    counter = 0
    while True:
        numerator = H - p*H(s)/p(s)  # H(s):Hλ(s), s:s_λ
        denominator = np.poly1d([1,-s])  # z-s
        H_new, remainder = np.polydiv(numerator,denominator)  # 多項式の除算
        H_new /= H_new.c[0]
           
        s += -p(s)/H_new(s)
        if abs(p(s)) <EPS or counter==LIMIT:
            break  
        
        H = H_new        
        counter += 1
    return s   

    
ang = plot_complex(Root_of_func)

以下に計算結果を示します。