DKA法

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

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

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

これを微分すると，

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

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

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

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

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

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

DKA法の手順†

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

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

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

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

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

/*!
 * Aberthの方法で初期値を算出
 * @param[out] z 解探索のための初期値
 * @param[in] c 多項式の係数
 * @param[in] n 方程式の次数
 * @param[in] max_iter 半径計算のための最大反復数
 * @param[in] eps 半径計算のための許容誤差
 * @return 1:成功,0:失敗
 */
int aberth(vector<complexf> &z, const vector<complexf> &c, int n, int max_iter, double eps)
{
    // 半径算出のための方程式の係数
    vector<complexf> cd(n+1, 1.0); // 係数c'
    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();
    // zの多項式をz+c1/nで割っていくことでwの多項式の係数を求める
    for(int i = n; i > 1; --i){
        horner(a, c1n, tmp, rm, i);
        cd[i] = rm;
        a = tmp;
    }
    cd[1] = a[1]+c1n;

    // 多項式S(w)の係数
    vector<double> b(cd.size());
    b[0] = abs(cd[0]);
    for(int i = 1; i <= n; ++i){
        b[i] = -abs(cd[i]);
    }

    // Aberthの初期値の半径をニュートン法で算出
    int m = 0; // 係数cの中で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);

    // Aberthの初期値
    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関数は以下．

/*!
  * ホーナー法(組立除法)
   - P(x) = a0 x^n + a1 x^(n-1) + ... + a_(n-1) x + a_n を (x-b) で割ったときの商と余りを返す
   - 商はn-1次の多項式の係数として返す
  * @param[in] a 代数方程式の係数
  * @param[in] b 割る1次式の係数(x-b)
  * @param[out] c 商であるn-1次の多項式の係数
  * @param[out] rm あまり
  * @param[in] n 方程式の次数(配列aの大きさはn+1,配列cはn)
  */
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];
    }
}

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

/*!
 * ワイヤストラス法(DK公式)
 * @param[inout] z 初期値位置を受け取り，解を返す
 * @param[in] c 多項式の係数
 * @param[in] n 方程式の次数
 * @param[inout] max_iter 最大反復数(反復終了後,実際の反復数を返す)
 * @param[inout] eps 許容誤差(反復終了後,実際の誤差を返す) 
 * @return 1:成功,0:失敗
 */
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;

        // DK式の計算
        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;
}

↑