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법을 사용해 보면 어떻게 될까요? 점화식의 우변 2항목의 분모는 $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. News.

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. 감의, W. 우안 d T. 야마모토, 리 MS 코큐 6의. 915 (1995) 225-249.