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

まず,代数方程式が複素解を持つとして,未知数をrf_dka.eq2.gifで表すことにする.

rf_dka.eq3.gif

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

rf_dka.eq5.gif

これを微分すると,

rf_dka.eq6.gif

rf_dka.eq7.gifが解rf_dka.eq8.gifに等しいとき,rf_dka.eq9.gifであるので,上式の左辺は第1項のみ残る. 同様にrf_dka.eq10.gifのとき左辺第rf_dka.eq11.gif項のみ残る.これを式で表すと,

rf_dka.eq12.gif

ここで,rf_dka.eq13.gifである. rf_dka.eq14.gifrf_dka.eq1.gifのそれぞれの解の近似値とすると,ニュートン法の公式から,

rf_dka.eq15.gif

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

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

rf_dka.eq17.gif

ここで,rf_dka.eq18.gifrf_dka.eq19.gifは複素平面上でrf_dka.eq1.gif個の解を包むような円の半径である. また,ここでのrf_dka.eq20.gifrf_dka.eq21.gifを示す.

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

DKA法の手順

  1. 中心rf_dka.eq22.gifと半径rf_dka.eq19.gifを計算し,初期値rf_dka.eq16.gifを計算
  2. DK式でrf_dka.eq14.gifを更新 収束するまで2を繰り返すことで,近似解を得る.

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

rf_dka.eq24.gif

ここで,rf_dka.eq25.gifはホーナーの方法で得られた係数である. rf_dka.eq26.gif でかつ,すべての係数rf_dka.eq25.gifが非ゼロとしたとき,rf_dka.eq27.gifの解は,

rf_dka.eq28.gif

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

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

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
/*!
 * 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> a;
    complexf zc = -c[1]/(c[0]*(double)n);
    horner(c, a, n, zc);
 
    vector<double> b(a.size());
    b[0] = abs(a[0]);
    for(int i = 1; i <= n; ++i){
        b[i] = -abs(a[i]);
    }
 
    // Aberthの初期値の半径をニュートン法で算出
    double r = 100.0;
    newton(b, n, r, max_iter, eps);
 
    // Aberthの初期値
    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;

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

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
/*!
 * ホーナー法(組立除法)の係数の値を求める
 * @param[in] a 代数方程式の係数
 * @param[out] b x=x0のときの組立除法の係数
 * @param[in] x0 係数を計算するxの値
 */
template<class T>
inline void horner(const vector<T> &a, vector<T> &b, int n, T x0)
{
    if(n <= 2) return;
 
    b = a;
    for(int i = 0; i <= n; ++i){
        for(int j = 1; j <= n-i; ++j){
            b[j] += x0*b[j-1];
        }
    }
}

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

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
/*!
 * ワイヤストラス法(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;
}

Aberthの初期値について

代数方程式の解と係数の関係を導いてみる. まず,2次方程式の場合を考える.

rf_dka.eq34.gif

rf_dka.eq35.gifを使うと方程式は以下のように書き直せる.

rf_dka.eq36.gif

元々の式と係数を比べることで,以下の関係式が得られる.

rf_dka.eq37.gif

同様にしてn次方程式の場合,

rf_dka.eq38.gif

よって,n次方程式の解と係数の関係式は,

rf_dka.eq39.gif

ここで,rf_dka.eq40.gifである. たとえば,3次方程式(rf_dka.eq41.gif)の場合,

rf_dka.eq42.gif

rf_dka.eq43.gifのときの関係式から,

rf_dka.eq44.gif

Aberthの初期値では解rf_dka.eq45.gifの複素平面における重心

rf_dka.eq46.gif

を中心とした半径rf_dka.eq19.gifの円周上に等間隔にrf_dka.eq16.gifを配置する.つまり,

rf_dka.eq47.gif

となる.オイラーの公式(rf_dka.eq48.gif)から,

rf_dka.eq49.gif

が得られる.ただし,ここでのrf_dka.eq50.gifrf_dka.eq21.gifを示す. また,rf_dka.eq51.gifを角度に加えたのは,初期値rf_dka.eq16.gifが実軸に対称に分布することを防ぐためである(共役複素数になることを防ぐ).


*1 http://en.wikipedia.org/wiki/Durand%E2%80%93Kerner_method
*2 O. Aberth, ``Iteration methods for finding all zeros of a polynomial simultaneously'', Math. Comput. 27, pp.339-344, 1973.

添付ファイル: filerf_dka.eq49.gif 639件 [詳細] filerf_dka.eq17.gif 617件 [詳細] filerf_dka.eq10.gif 666件 [詳細] filerf_dka.eq9.gif 588件 [詳細] filerf_dka.eq8.gif 606件 [詳細] filerf_dka.eq7.gif 635件 [詳細] filerf_dka.eq15.gif 604件 [詳細] filerf_dka.eq21.gif 604件 [詳細] filerf_dka.eq29.gif 525件 [詳細] filerf_dka.eq37.gif 610件 [詳細] filerf_dka.eq1.gif 575件 [詳細] filerf_dka.eq33.gif 553件 [詳細] filerf_dka.eq20.gif 571件 [詳細] filerf_dka.eq6.gif 641件 [詳細] filerf_dka.eq45.gif 640件 [詳細] filerf_dka.eq14.gif 604件 [詳細] filerf_dka.eq48.gif 556件 [詳細] filerf_dka.eq38.gif 584件 [詳細] filerf_dka.eq23.gif 649件 [詳細] filerf_dka.eq30.gif 635件 [詳細] filerf_dka.eq35.gif 559件 [詳細] filerf_dka.eq11.gif 604件 [詳細] filerf_dka.eq32.gif 597件 [詳細] filerf_dka.eq47.gif 616件 [詳細] filerf_dka.eq25.gif 645件 [詳細] filerf_dka.eq18.gif 672件 [詳細] filerf_dka.eq16.gif 477件 [詳細] filerf_dka.eq27.gif 602件 [詳細] filerf_dka.eq24.gif 564件 [詳細] filerf_dka.eq44.gif 553件 [詳細] filerf_dka.eq22.gif 559件 [詳細] filerf_dka.eq39.gif 608件 [詳細] filerf_dka.eq13.gif 629件 [詳細] filerf_dka.eq34.gif 642件 [詳細] filerf_dka.eq40.gif 561件 [詳細] filerf_dka.eq42.gif 600件 [詳細] filerf_dka.eq36.gif 529件 [詳細] filerf_dka.eq4.gif 545件 [詳細] filerf_dka.eq28.gif 601件 [詳細] filerf_dka.eq12.gif 641件 [詳細] filerf_dka.eq51.gif 603件 [詳細] filerf_dka.eq19.gif 605件 [詳細] filerf_dka.eq2.gif 667件 [詳細] filerf_dka.eq43.gif 471件 [詳細] filerf_dka.eq26.gif 623件 [詳細] filerf_dka.eq5.gif 644件 [詳細] filerf_dka.eq3.gif 665件 [詳細] filerf_dka.eq31.gif 591件 [詳細] filerf_dka.eq46.gif 586件 [詳細] filerf_dka.eq50.gif 651件 [詳細] filerf_dka.eq41.gif 673件 [詳細]

トップ   編集 凍結 差分 バックアップ 添付 複製 名前変更 リロード   新規 一覧 単語検索 最終更新   ヘルプ   最終更新のRSS
Last-modified: 2012-06-27 (水) 12:15:32 (2641d)