DKA法

代数方程式は重根も別々に考えると個の複素数解を持つ(ここでは実数も複素数の一部としている)．ここまでの方法では適切な初期値や範囲を与えることで1つの解が得られた．一方で代数方程式での複数解を求める場合は初期値や探索範囲を適切に設定し直さなければならない．ここでは個の複素数解を一度に計算する方法について述べる．

まず，代数方程式が複素解を持つとして，未知数をで表すことにする．

代数方程式におけるニュートン法を考える．上記の代数方程式の解をとすると，代数方程式は以下のように変形できる．

これを微分すると，

が解に等しいとき，であるので，上式の左辺は第1項のみ残る．同様にのとき左辺第項のみ残る．これを式で表すと，

ここで，である．をのそれぞれの解の近似値とすると，ニュートン法の公式から，

ここで，．この式をワイヤストラス(Weierstrass)法，もしくは，デュラン・ケルナー(Durand-Kerner)法と呼ぶ *1．以下では上式をDK式と呼ぶことにする．

DK式により個の複素数解が計算できるが，ニュートン法を用いているための初期値を適切に選ぶ必要がある．そこで，以下のアバース(Aberth)の初期値を用いる．

ここで，，は複素平面上で個の解を包むような円の半径である．また，ここでのはを示す．

Aberthの初期値を用い，DK式を反復計算することで，個の根を得る方法をDKA(Durand-Kerner-Aberth)法と呼ぶ．

DKA法の手順†

中心と半径を計算し，初期値を計算
DK式でを更新収束するまで2を繰り返すことで，近似解を得る．

Aberthの初期値の半径を求める方法は様々に提案されているが，ここでは Aberthの方法*2を説明する．として以下の方程式を得る．

ここで，はホーナーの方法を応用した組立除法で得られた係数である．でかつ，すべての係数が非ゼロとしたとき，の解は，

の代数方程式を解いた時の解を半径とした円内に存在する．よって，の解は，内にある．このを初期値を算出するための半径として用いる．は単一の実数解を持つのでニュートン法を使って解ける．

DKA法を使って複数の根を一度に計算するコード例を以下に示す．まず，Aberthの初期値を求める関数である．


int aberth(vector<complexf> &z, const vector<complexf> &c, int n, int max_iter, double eps)
{
        vector<complexf> cd(n+1, 1.0);     double c1n = -c[1].real()/n; 
    cd[0] = c[0];
    double rm;     vector<double> a(n+1), tmp(n);     for(int i = 0; i <= n; ++i) a[i] = c[i].real();
        for(int i = n; i > 1; --i){
        horner(a, c1n, tmp, rm, i);
        cd[i] = rm;
        a = tmp;
    }
    cd[1] = a[1]+c1n;
 
        vector<double> b(cd.size());
    b[0] = abs(cd[0]);
    for(int i = 1; i <= n; ++i){
        b[i] = -abs(cd[i]);
    }
 
        int m = 0;     for(int i = 0; i <= n; ++i) m += (fabs(c[i].real()) > 1e-6 ? 1 : 0);
    double rmax = 0.0;     for(int i = 1; i <= n; ++i){
        double ri = pow(m*fabs(c[i].real())/fabs(c[0].real()), 1.0/(i+1.0));
        if(ri > rmax) rmax = ri;
    }
    double r = rmax;
    newton(b, n, r, max_iter, eps);
 
        complexf zc = -c[1]/(c[0]*(double)n);
    for(int j = 0; j < n; ++j){
        double theta = (2*RX_PI/(double)n)*j+RX_PI/(2.0*n);
        z[j] = zc+r*complexf(cos(theta), sin(theta));
    }
 
    return 1;
}

ここでは複素数を扱うために，C++のcomplexを使っている．

#include <complex>	// 複素数
typedef complex<double> complexf;

半径を算出するために，ニュートン法のコードを変更し，係数列を与えると関数値と導関数値をホーナー法で計算するようにした．また，半径算出のためのの係数算出に使っているhorner関数は以下．

 
inline void horner(const vector<double> &a, double b, vector<double> &c, double &rm, int n)
{
    if(n <= 1) return;
    rm = a[0];     c.resize(n);
    for(int i = 1; i < n+1; ++i){
        c[i-1] = rm;
        rm *= b;
        rm += a[i];
    }
}

初期値を算出後，以下のワイヤストラス法の関数に多項式の係数とともに渡すことで解を得る．


int weierstrass(vector<complexf> &z, const vector<complexf> &c, int n, int &max_iter, double &eps)
{
    double e = 0.0, ej;
 
    vector<complexf> zp;
    complexf f, df;
    int k;
    for(k = 0; k < max_iter; ++k){
        zp = z;
 
                for(int j = 0; j < n; ++j){
            f = func(z[j], c, n);
            df = c[0];
            for(int i = 0; i < n; ++i){
                if(i != j){
                    df *= zp[j]-zp[i];
                }
            }
 
            z[j] = zp[j]-f/df;
        }
 
                e = 0.0;
        for(int j = 0; j < n; ++j){
            if((ej = abs(func(z[j], c, n))) > e){
                e = ej;
            }
        }
 
                if(e < eps){
            max_iter = k; eps = e;
            return 1;
        }
    }
 
    eps = e;
    return 0;
}