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

1. ラゲール法(Laguerre’s Method)の導出
2. ラゲール法の計算プログラム(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

Laguerre’s method – Wikipedia

関連記事：

ニュートン法による複素求根アルゴリズム~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)

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