Durand-Kerner-Aberth法
$a_1, a_2, \cdots, a_n$を複素数として、以下のような多項式
P(z) = z^n + a_1 z^{n-1} + \cdots + a_{n-1}z + a_n
の零点を探索する問題を考えます。Durand-Kerner-Aberth (DKA)法は、すべての零点を同時に求めることができるエレガントなアルゴリズムです。DKA法はある非線形方程式に対するNewton法と等価であることが知られており、零点に縮退がなければ、少なくとも局所的には収束します。原型はWeierstrassが1903年に既に発見していたらしいです。
アルゴリズム
(1) $z_1^0, z_2^0, \cdots, z_n^0$ を以下の値にセットします。これはAberthの初期値と呼ばれます。
z_j^0 = -\frac{a_1}{n} + r e^{\pi i\theta_j/n}, \quad
\theta_j = 2(j-1)+\frac{1}{2}
ここで、$r$ は適当な実数です。
(2) $z_j^1, z_j^2, \cdots$ を以下の漸化式
z_j^{l+1} = z_j^l - \frac{P(z_j^l)}{\prod_{j\neq k} (z_j^l - z_k^l)}, \quad l\geq 0
にしたがって計算します。適当な回数反復すると、$z_j^l$ は $P(z)$ の零点に収束します。
Juliaでの実装例
上の式を愚直に書いていきます。
まず、Aberthの初期値を与えるコードを示します。$z_1^0, \cdots, z_n^0 $ のリストが戻り値になっています。
function initial_condition(a₁, n, r)
θⱼ = [2(j-1)+1/2 for j in 1:n]
return r*exp.(1.0im*π*θⱼ/n) .- a₁/n
end
ここで、exp.(a) はelement-wiseに演算することを表します。例えば、
exp.([a₁, a₂, a₃]) == [exp(a₁), exp(a₂), exp(a₃)]
はtrueになります。exp() に限らず任意の関数について、同様のことができます。ちょっと分かりにくいですが、ベクトルとスカラー量の演算に対して.を使った場合は
[a₁, a₂, a₃] .- b₁ == [a₁-b₁, a₂-b₁, a₃-b₁]
のようになります。
次に漸化式を書きます。
# a = [a₁, a₂, ...]
# num: number of iterations
function iteration(a, z, num)
n = length(a)
P(z) = sum(a .* [z^k for k in n-1:-1:0]) + z^n
for k in 1:num
denom = [prod([z[i]-z[j] for j in 1:n if i != j]) for i in 1:n]
z = z .- P.(z) ./ denom
end
return z
end
最後に、上の二つをまとめます。多項式の係数、初期値のコントロールパラメータ、反復回数を与えると、零点の候補が返ってきます。
function find_roots(a, r, num)
n = length(a)
z = initial_condition(a[1], n, r)
return iteration(a, z, num)
end
以上で完成です。さっそく試して見ましょう。
（例1）$P(z) = (z-2)(z+4) = z^2 + 2z - 8$
$r = 10$, すなわち、厳密解の外側に初期値をセットします。
julia> find_roots([2, -8], 10, 10) ≈ [2, -4]
true
次の図は、複素平面上で $z_j^l$ がどのように動いているかを示したものです。星印は厳密な零点を表します。
今度は $r = 1$, すなわち、厳密解の内側に初期値をセットします。
julia> find_roots([2, -8], 1, 10) ≈ [2, -4]
true
（例2）
P(z)
= (z-2)(z+4)(z-i)(z+3i) \\
= z^4 + (2 + 2i) z^3 + (-5 + 4i) z^2 + (6 - 16i) z - 24
julia> find_roots([2+2im, -5+4im, 6-16im, -24], 10, 10) ≈ [2, 1im, -4, -3im]
true
縮退がある場合（実験）
零点に縮退がある場合にあえてDKA法を使ってみるとどうなるでしょうか？漸化式の右辺二項目の分母は $z_j^l - z_k^l$ という形をしていることから、$z_j^l$ たちがある程度互いに接近すると、反発することが予想されます。この様子を観察してみましょう。
(例3) $P(z) = (z-1)^{10}$
まとめ
DKA法を紹介しました。アニメーションの描画を含むサンプルコードは、こちらで公開しています。
参考文献
-
Weierstrass, Neuer Beweis des Satzes, dass jede ganze rationale Funktion einer Veraenderderlichen dargestellt werden kann als ein Product aus lineare Funktionen derselben Veraenderlichen. Gesammelte Werke 3 (1903), Johnson New York, 1967, 251-269.
-
I. E. Durand, Solutions Numerique des Equations Algebriques. Tome I: Equations du Type $F(x)=0$. Racines d’une Polynome, Masson, Paris, 1960, 279-281.
-
I. O. Kerner, Ein Gesamtschrittverfahren zur Berechnung der Nullstelle von Polynomen. Numer. Math. 8(1966), 290-294.
-
O. Aberth, Iteration methods for finding all zeros of a polynomial simultaneously. Math. Comp. 27(1973), 339-344.
-
S. Kanno, W. Liu and T. Yamamoto, RIMS Kokyuroku No.915 (1995) 225-249.