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

ラゲール法(Laguerre’s Method)の導出
求根アルゴリズムのラゲール法(Laguerre’s Method)とは、\(n\)次の多項式\(f(z)\)の根、すなわち\(f(z)=0\)の解\(z\)を求める手法のひとつである。
\(f(z)=0\)が解、
$$z = z_1, z_2, z_3, \cdots, z_n $$
を持つとすると、以下のように書ける。(\(C\)は係数)
$$ f(z) = C (z-z_1) (z-z_2) (z-z_3) \cdots (z-z_n) $$
ここで、以下の関数\(G(z),H(z)\)を考える。
$$ \begin{eqnarray}
G(z) &=& \frac{d}{dz} \ln{f(z)} \\
&=& \frac{f'(z)}{f(z)} \tag{1-1} \\
&=& \sum_{i=1}^n \frac{1}{z-z_i} \tag{1-2} \\
\\
H(z) &=& – \frac{dG(z)}{dz} \\
&=& \frac{ (f'(z))^2-f(z)f”(z) }{ (f(z))^2 } \tag{2-1} \\
&=& \sum_{i=1}^n \frac{1}{(z-z_i)^2} \tag{2-2}
\end{eqnarray} $$
これから求めていく\(f(z)\)の根を\(z_1\)と設定して以下話を進める。
\( \alpha(z), \beta(z), \delta_i(z) \)を下記のようにおいてみる。
$$
\alpha (z) = \frac {1}{z-z_1} \tag{3} \\
\beta(z) + \delta_i(z) = \frac{1}{z-z_i} \ \ \ (i = 2,3,\cdots , n)
$$
ただし\(\beta(z) \)は \( 1/(z-z_i) (i=2,3,\cdots,n)\)の平均値である。
$$ \beta(z)=\frac{1}{n-1}\sum_{i=2}^n \frac{1}{z-z_i} $$
これより、明らかに\( \delta_i (z)\) の \( i=2,3,\cdots,n \)における和は\(0\)である。
$$ \sum_{i=2}^n \delta_i(z) = 0 $$
さらに、以下のように\( \Delta^2(z) \)を導入する。
$$ \Delta^2(z) = \sum_{i=2}^n \delta_i^2 (z) $$
ここで、\( \alpha(z), \beta(z), \Delta^2(z) \)を用いて式(1-2), (2-2)の\( G(z), H(z) \)を表してみる。
$$ G(z) = \alpha(z) + (n-1) \beta(z) $$
$$ H(z) = \alpha^2(z) + (n-1) \beta^2(z) + \Delta^2(z) $$
上式から、\( \alpha(z) \)を求める。
$$ \alpha(z) = \frac{ G(z) \pm \sqrt{ (n-1)[nH(z) – G^2(z) -n \Delta^2(z)] } }{ n } $$
式(3)とより
$$ z_1 = z \ – \frac{n}{ G(z) \pm \sqrt{ (n-1)[nH(z) – G^2(z) -n \Delta^2(z)] } } $$
上式は\( z_1 \)についての正確な表現である。\(z\)が\(z_1\)に近い場合には\( \Delta^2(z) \)は他の項と比較して小さくなるので、その項を無視する。
そうするとラゲール法のアルゴリズムが得られる。
$$ z_1^{(k+1)} = z_1^{(k)} – \frac{n}{ G(z_1^{(k)}) \pm \sqrt{ (n-1)[nH(z_1^{(k)}) – G^2(z_1^{(k)}) ] } } \tag{4}$$
分母の符号は、絶対値が大きくなるほうを採用する。
初期値\(z_1^{(0)}\)を式(4)に代入して逐次的に\(z_1^{(k)} \ \ (k=1,2,3\cdots)\)を求めていき、収束するまでこれを繰り返す。
\(H(z_1^{(k)}), G(z_1^{(k)})\)を計算するには式(1-1),(2-1)を使用する。
参考:
E. HANSEN et al., “SOME MODIFICATIONS OF LAGUERRE’S METHOD” BIT 17 (1977) 409-417
関連記事:
ニュートン法による複素求根アルゴリズム~Pythonコード付き~
Jenkins-Traubアルゴリズムで多項式の根を求める~Pythonコード付き~
ラゲール法の計算プログラム(Python)
以下のプログラムでは、
$$f(z) = z^5 +1 $$
を例に、初期値\(z_0\)を\(\{ z \mid -10<Re(z)<10, -10<Im(z)<10 \} \)の範囲でラゲール法を適用し、根の偏角の大きさをカラーマップで表示しています。
# -*- coding: utf-8 -*- """ 複素平面上の値を初期値として、方程式func(z)=0の根zをラゲール法で計算し、 その偏角の大きさを2次元カラーマップする。 """ import numpy as np import matplotlib.pyplot as plt LIMIT = 100 # 計算回数の上限 EPS = 1E-10 # ゼロ点の判定範囲. abs(f(z)) < eps ならばゼロ点とする. # 根を求めたい多項式 f = np.poly1d([1,0,0,0,0,1]) # f(z)=z^5+1 df = f.deriv() # 1階微分 ddf = f.deriv(2) # 2階微分 n = f.order # 多項式f(z)の次数 print(f) # 初期値の計算範囲 real_min=-10 real_max=-real_min im_min=-10 im_max=-im_min N_real=501 # ステップ数 N_im=501 # ステップ数 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) plt.xlabel("Real part") plt.ylabel("Imaginary part") 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 def Root_of_func(z): # ラゲール法アルゴリズム counter = 0 while True: if abs(f(z)) < EPS: break G = df(z)/f(z) H = G**2 - ddf(z)/f(z) aDenom_pls = G + np.sqrt((n - 1)*(n*H - G**2)) # プラス aDenom_mns = G - np.sqrt((n - 1)*(n*H - G**2)) # マイナス if aDenom_pls*aDenom_mns == 0 : break # 絶対値が大きいほうを採用 a = n/aDenom_pls if abs(aDenom_pls)>=abs(aDenom_mns) else n/aDenom_mns z -= a # 収束判定 counter += 1 if abs(a) < EPS: #abs(real_max - real_min)*1E-10: break if counter == LIMIT: break return z ang = plot_complex(Root_of_func)
以下に計算結果を示します。

ニュートンフラクタルの場合とは異なる図形が現れました。その他の方程式でも様々な図形が現れて非常に面白いので、上記のPythonプログラムをコピペして、いろいろ試してみてください。
-
前の記事
これが解けたらIQ150!?~囚人の帽子の問題~【TED-Ed】 2021.03.14
-
次の記事
Jenkins-Traubアルゴリズムで多項式の根を求める~Pythonコード付き~ 2021.03.20