WANtaroHP (DKA)

DKA 法

概要

高次の代数方程式の重解を含む全ての解を求めるための方法として，DKA 法を紹介する．

これは，Newton 法を代数方程式の解法に発展させた Durand-Kerner の公式を用い，Newton 法適用のための解の初期値を与えるため Aberth の初期値を用いる方法である．重解の場合の収束性は良くない (4倍精度計算にしないと使い物にならない) が，複素数解を含む全ての解を一発で求めることができ，便利である．

DK 法の理論

以下の方程式を考える．

\begin{equation} f(z)=a_0 z^n+a_1 z^{n-1}+\cdots+a_{n-1} z+a_n=0 \end{equation}

上式のn個の真の解を $\alpha_j (j=1\sim n)$ とすれば，

\begin{equation} f(z)=a_0 (z-\alpha_1)(z-\alpha_2)\cdots(z-\alpha_n)=a_0 \prod_{\substack{j=1}}^n (z-\alpha_j) \end{equation}

これを微分すれば，

\begin{equation} f'(z)=a_0 \prod_{\substack{j=2}}^n (z-\alpha_j) +a_0 \prod_{\substack{j=1\\j\neq 2}}^n (z-\alpha_j)+\cdots +a_0 \prod_{\substack{j=1}}^{n-1} (z-\alpha_j) \end{equation}

ここで，例えば， $z=\alpha_1$ だとすれば，上式右辺第1項を除いて $\alpha_1-\alpha_1=0$ となる項を含むため，

\begin{equation} f'(\alpha_1)=a_0 \prod_{\substack{j=2}}^n (\alpha_1-\alpha_j) \end{equation}

また， $z=\alpha_2$ だとすれば，上式右辺第2項を除いて $\alpha_2-\alpha_2=0$ となる項を含むため，

\begin{equation} f'(\alpha_2)=a_0 \prod_{\substack{j=1\\j\neq 2}}^n (\alpha_2-\alpha_j) \end{equation}

となり，結局，1項のみしか残らず，微分した値は，

\begin{equation} f'(\alpha_i)=a_0 \prod_{\substack{j=1\\j\neq i}}^n (\alpha_i-\alpha_j) \end{equation}

となる．

ここで， $\alpha_i$ の近似値を $z_i$ とし， $f'(z_i)$ を以下のとおり近似する．

\begin{equation} f'(z_i)\doteqdot a_0 \prod_{\substack{j=1\\j\neq i}}^n (z_i-z_j) \end{equation}

以上より，Newton 法の考え方を適用し，以下に示す Durand-Kerner の公式を導くことができる．

\begin{equation} z_i^{(k+1)}=z_i^{(k)}-\cfrac{f(z_i^{(k)})}{a_0 \prod_{\substack{j=1\\j\neq i}}^n (z_i^{(k)}-z_j^{(k)})} \end{equation}

ここに $z_i$ は，方程式 $f(z)=0$ の $i$ 番目の解の近似値である．

Aberth の初期値

また，Newton 法を適用するために必要となる初期値は，以下に示す，Aberth の初期値を用いることができる．

\begin{equation} z_j^{(0)}=-\cfrac{a_1}{n a_0}+R\cdot\exp\left\{i\cdot\cfrac{2\pi}{n}\left(j-\cfrac{3}{4}\right)\right\} \qquad i=\sqrt{-1} \end{equation}

上式は，代数方程式の解と係数の関係

\begin{equation} \cfrac{\alpha_1+\alpha_2+\cdots+\alpha_n}{n}=-\cfrac{a_1}{n a_0} \end{equation}

より，複素平面において，解の重心 $(-a_1/(n a_0), 0)$ を中心とした半径 $R$ の円の中に解が存在するよう $R$ を設定し，その円弧上に均等に初期値を配置しようとするものである．

方程式の書き換え

ここで解きたい方程式を以下のように書き換える．

\begin{equation} g(z)=z^n+c_1 z^{n-1}+\cdots+c_{n-1} z+c_n=0 \qquad c_i=a_i/a_0 \end{equation}

上式より Durand-Kerner の公式，Aberth の初期値は以下のとおり現される．

\begin{equation} z_i^{(k+1)}=z_i^{(k)}-\cfrac{f(z_i^{(k)})}{\prod_{\substack{j=1\\j\neq i}}^n (z_i^{(k)}-z_j^{(k)})} \end{equation}

\begin{equation} z_j^{(0)}=-\cfrac{c_1}{n}+R\cdot\exp\left\{i\cdot\cfrac{2\pi}{n}\left(j-\cfrac{3}{4}\right)\right\} \qquad i=\sqrt{-1} \end{equation}

Aberth の初期値での $R$ は，以下の方程式の正の実数解を用いることができる．( $n-1$ 項が無い事に注意)

\begin{equation} x^n-|c_2| x^{n-2}-|c_3| x^{n-3}-\cdots-|c_{n-1}| x-|c_n|=0 \end{equation}

$c_n\neq 0$ とすれば，上式は $x=0$ で負の値をとるため， $x$ が正の領域で実数解が存在し，これは Newton 法により求めることができる．

参考

Newton 法

上記過程は，Aberth の初期値の $R$ の設定を除いて，全て複素数計算で行う必要があるが，ここでは単純な事例として，代数方程式の1つの実数解を得るための標準的な Newton 法の手順を示しておく．

代数方程式の1つの実数解を得るための Newton 法の標準手順

\begin{equation*} f(x) = a_0 x^n + a_1 x^{n-1}+\cdots+ a_{n-2} x^2+a_{n-1} x + a_n \end{equation*}

\begin{equation*} f'(x)=n a_0 x^{n-1}+(n-1) a_1 x^{n-2}+\cdots+2 a_{n-2} x +a_{n-1} \end{equation*}

`xx`	変数`x`	`x0`	変数の初期値
`a(i)`	代数方程式係数`a_i`	`ii`	収束計算用カウンタ
`ff`	関数値`f(x)`	`iimax`	収束計算最大繰返数
`df`	関数微分値`f'(x)`	`eps`	収束判定用の絶対値が小さな数値

! main routine
xx=x0
do ii=1,iimax
    ff=a(0)
    do i=1,n
        ff=xx*ff+a(i)
    end do
    df=n*a(0)
    do i=1,n-1
        df=xx*df+(n-i)*a(i)
    end do
    xx=xx-ff/df
    if(abs(ff/df) < eps)exit
end do

複素数の演算

プログラムにおいて複素数タイプを用いずに複素数の演算を行うには，複素数を実数部と虚数部に分けて演算することが必要になる．そのための基本的な演算公式を示す．

複素数 $f$, $z$, $p$ を実数部と虚数部に分けて以下のように定義する．，

\begin{equation*} f=f_1+f_2\cdot i \qquad z=z_1+z_2\cdot i \qquad p=p_1+p_2\cdot i \end{equation*}

複素数の乗算および除算は以下のとおり表現できる．

\begin{align*} &f*z &=&(f_1 z_1-f_2 z_2)+(f_2 z_1+f_1 z_2)\cdot i \\ & \\ &p*(z(i)-z(j))&=&p_1(z_1(i)-z_1(j))-p_2(z_2(i)-z_2(j)) \\ & &~&+\{p_1(z_2(i)-z_2(j))+p_2(z_1(i)-z_1(j))\} \\ & \\ &\cfrac{f}{p} &=&\cfrac{f_1 p_1+f_2 p_2}{p_1^2+p_2^2}+\cfrac{f_2 p_1-f_1 p_2}{p_1^2+p_2^2}\cdot i \end{align*}

Programs

Filename	Description
f90_DKA.txt	program for solution of algebric equations by DKA method
inp_test10.txt	input data sample (1) equation of the 10th degree
inp_test30.txt	input data sample (2) equation of the 30th degree
inp_test70.txt	input data sample (3) equation of the 7th degree
out_test10.txt	output data sample (1) equation of the 10th degree
out_test30.txt	output data sample (2) equation of the 30th degree
out_test70.txt	output data sample (3) equation of the 7th degree