非線形方程式の根を数値的に求める方法について



線形方程式と非線型方程式

rootfinding.eq1.gifについての方程式rootfinding.eq2.gifの解を求めたい. もし,rootfinding.eq3.gif

rootfinding.eq4.gif

のような線形方程式(1次方程式)であるならば,その解は,

rootfinding.eq5.gif

と計算できる(ただし,rootfinding.eq6.gif). しかし,rootfinding.eq3.gifが非線形方程式であったならどうだろうか. たとえば,以下のrootfinding.eq7.gif次代数方程式(algebraic equation)を考える.

rootfinding.eq8.gif

2次方程式(rootfinding.eq9.gif)ならば解の方程式,3次や4次方程式の場合はカルダノ(Cardano)やフェラリ(Ferrari)の公式で 解析的に解くことができるが,5次以上だと直接的な解法はない.

また,超越方程式(transcendental equation)と呼ばれる無限次代数方程式(rootfinding.eq10.gif)もある. たとえば,rootfinding.eq11.gifは無限べき級数で表されるため,rootfinding.eq11.gifを含む方程式は超越方程式となる. 超越方程式に関してもほとんどの場合直接的な解法がない.

以下では非線形方程式を数値的に解く方法として,

  • 2分法
  • ニュートン・ラフソン法
  • 多次元のニュートン・ラフソン法
  • ホーナー法
  • DKA法

について述べる.

2分法

関数rf_bisection.eq1.gifを区間rf_bisection.eq2.gifを考えたとき, rf_bisection.eq3.gifrf_bisection.eq4.gifの符号が異なる場合(rf_bisection.eq5.gif),区間rf_bisection.eq2.gif内に少なくとも1つの解が存在する. 区間rf_bisection.eq2.gifの中点での関数値rf_bisection.eq6.gifを求め,その値の符号がrf_bisection.eq4.gifと同じならば, 解はrf_bisection.eq7.gif内に存在する. このとき,解の存在する区間はrf_bisection.eq8.gifになっている. 同様にして,rf_bisection.eq7.gifの中点での関数値を使うことで,解の存在区間はrf_bisection.eq9.gifとなる. この作業を繰り返すことで方程式の解を求める方法を2分法(bisection method)と呼ぶ.

bisection.jpg
図1 2分法

図1は2分法で近似解を求める手順を示している. 現在の解の存在区間をrf_bisection.eq10.gifとする. rf_bisection.eq10.gifrf_bisection.eq2.gifで初期化される. 中点をrf_bisection.eq11.gifとする. 図1ではrf_bisection.eq12.gifであるので,解はrf_bisection.eq13.gifに存在するので, rf_bisection.eq14.gifとする. 次の中点をrf_bisection.eq15.gifであり, 図1ではrf_bisection.eq16.gifであるので,解はrf_bisection.eq17.gifに存在する. これを反復して,rf_bisection.eq18.gifと求めていく.

2分法の手順を以下に示す.

  1. rf_bisection.eq19.gif, rf_bisection.eq20.gifとする.
  2. 中点rf_bisection.eq21.gifにおける関数値rf_bisection.eq22.gifを求める.
  3. rf_bisection.eq23.gifの場合rf_bisection.eq24.gifrf_bisection.eq25.gifの場合rf_bisection.eq26.gifとする.
  4. 2に戻る

反復を終了するための収束条件は以下となる.

rf_bisection.eq27.gif

2分法のコード例を以下に示す.

  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
/*!
 * 2分法(bisection method)
 * @param[in] func 関数値を与える関数ポインタ
 * @param[in] xl,xr 探索範囲
 * @param[out] x 解
 * @param[inout] max_iter 最大反復数(反復終了後,実際の反復数を返す)
 * @param[inout] eps 許容誤差(反復終了後,実際の誤差を返す) 
 * @return 1:成功,0:失敗
 */
int bisection(double func(const double), double xl, double xr, double &x, int &max_iter, double &eps)
{
    double f = func(xl);
    double fmid = func(xr);
 
    // 探索範囲の境界での値の符号が異なる場合のみ探索
    if(f*fmid >= 0.0) return 0.0;
 
    int k;
    double dx = fabs(xr-xl), xmid;
    for(k = 0; k < max_iter; ++k){
        xmid = 0.5*(xl+xr);    // 中点
        dx *= 0.5;
 
        // 中点での関数値を求める
        fmid = func(xmid);
 
        // 収束判定
        if(dx < eps || fmid == 0.0){
            x = xmid;
            max_iter = k; eps = dx;
            return 1;
        }
 
        // 新しい区間
        if(f*fmid < 0){
            xr = xmid;
        }
        else{
            xl = xmid;
            f = fmid;
        }
    }
 
    max_iter = k; eps = dx;
    return 0;
}

ニュートン・ラフソン法

2分法は指定した区間内の解は必ず求めるものの,解の収束は遅い. 関数rf_newton.eq1.gifが連続で微分可能関数であるならば,関数の勾配値を用いることで収束を早めることができる.

解の初期値をrf_newton.eq2.gifとする. rf_newton.eq3.gifでの関数rf_newton.eq1.gifの接線は,

rf_newton.eq4.gif

となる.この接線がrf_newton.eq5.gif軸と交差する(rf_newton.eq6.gif)点は,

rf_newton.eq7.gif

rf_newton.eq8.gifを次の近似値として同様の処理を繰り返すことで解に近づいていく.図2はこの手順を示したものである. この方法をニュートン・ラフソン法(Newton-Raphson method),もしくは単にニュートン法(Newton method)と呼ぶ. ニュートン法の式は以下となる.

rf_newton.eq9.gif

反復回数を多くすると限りなくrf_newton.eq10.gifの解に近づく. 反復を終了するための収束条件は,

rf_newton.eq11.gif

となる.

ニュートン法は2分法と異なり解のある区間でなく初期近似値1つだけを指定するだけでよく,収束も速い. しかし,rf_newton.eq1.gifだけでなくその微分値rf_newton.eq12.gifも必要で,初期値によっては解へ収束しない場合もある.

newton.jpg
図2 ニュートン法

ニュートン法のコード例を以下に示す.

  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
/*!
 * ニュートン・ラフソン法(Newton-Raphson method)
 * @param[in] func 関数値を与える関数ポインタ
 * @param[inout] x 探索開始位置を受け取り,解を返す
 * @param[inout] max_iter 最大反復数(反復終了後,実際の反復数を返す)
 * @param[inout] eps 許容誤差(反復終了後,実際の誤差を返す) 
 * @return 1:成功,0:失敗
 */
int newton(double func(const double), double dfunc(const double), 
           double &x, int &max_iter, double &eps)
{
    double f, df, dx;
    int k;
    for(k = 0; k < max_iter; ++k){
        f = func(x);
        df = dfunc(x);
        x = xn-f/df;
 
        // 収束判定
        dx = fabs(f/df);
        if(dx < eps || fabs(f) < eps){
            max_iter = k; eps = dx;
            return 1;
        }
    }
 
    max_iter = k; eps = dx;
    return 0;
}

多次元のニュートン・ラフソン法

ニュートン法を連立非線形方程式に一般化する.

rf_newton_md.eq1.gif

ベクトル表記では,

rf_newton_md.eq2.gif

ここで,

rf_newton_md.eq3.gif

である.

rf_newton_md.eq4.gifステップ目の近似値をrf_newton_md.eq5.gifとし, rf_newton_md.eq5.gifの周りで上式をテイラー展開する.

rf_newton_md.eq6.gif

ここで,rf_newton_md.eq7.gifrf_newton_md.eq8.gifrf_newton_md.eq9.gifの要素とするヤコビ行列である. 2次以上の項を無視すると,連立非線形方程式は以下となる.

rf_newton_md.eq10.gif

rf_newton_md.eq11.gifとすると,

rf_newton_md.eq12.gif

が得られる. この式はrf_newton_md.eq13.gifを未知数とした線形連立方程式であり, LU分解などで解くことで,rf_newton_md.eq13.gifが得られる. そして,以下の式でrf_newton_md.eq14.gifを計算する.

rf_newton_md.eq15.gif

例)2元連立非線形方程式の場合

rf_newton_md.eq16.gif

そして,

rf_newton_md.eq17.gif

である.よって,rf_newton_md.eq18.gifに関する式は以下となる.

rf_newton_md.eq19.gif

rf_newton_md.eq13.gifについて解くと,

rf_newton_md.eq20.gif

これらを用いてrf_newton_md.eq21.gifを更新する.

rf_newton_md.eq22.gif

ホーナー法(組立除法)

2分法やニュートン法を使って以下のrf_horner.eq1.gif次の代数方程式を解くときの演算回数を考える.

rf_horner.eq2.gif

rf_horner.eq3.gifの値を計算する際の乗算回数はrf_horner.eq4.gif回で 加算がrf_horner.eq1.gif回である. 乗算回数がrf_horner.eq1.gifの2乗で増えるため,たとえば,rf_horner.eq5.gifでは55回だが,rf_horner.eq6.gifだと5050回にもなる.

乗算回数を少なくするために代数方程式を以下のように変形する.

rf_horner.eq7.gif

この場合,一番内側の括弧内から計算していくと,

rf_horner.eq8.gif

となり,最終的にrf_horner.eq9.gifとなる. このときの乗算回数はrf_horner.eq1.gif回,加算もrf_horner.eq1.gif回である. 特にrf_horner.eq1.gifが大きな時に計算回数を大幅に減らすことができる. たとえば,乗算回数はrf_horner.eq5.gifで10回,rf_horner.eq6.gifでも100回で済む. この計算方法はホーナー(Horner)法と呼ばれる.

ホーナー法で代数方程式の値とその導関数を求めるコード例を以下に示す.

  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
/*!
 * ホーナー法で代数方程式の値を計算
 * @param[in] x 変数
 * @param[in] b 係数
 * @param[in] n 方程式の次数
 * @return 代数方程式の値
 */
template<class T>
inline T func_h(double x, const vector<T> &b, int n)
{
    T f = b[0];
    for(int i = 1; i <= n; ++i){
        f = b[i]+f*x;
    }
    return f;
}
/*!
 * ホーナー法で代数方程式の導関数値を計算
 * @param[in] x 変数
 * @param[in] b 係数
 * @param[in] n 方程式の次数
 * @return 代数方程式の導関数値
 */
template<class T>
inline T dfunc_h(double x, const vector<T> &b, int n)
{
    T df = n*b[0];
    for(int i = 1; i <= n-1; ++i){
        df = (n-i)*b[i]+df*x;
    }
    return df;
}

テンプレート関数にしているのはDKA法で複素数を扱う関係上,様々な型に対応できるようにしたいためである.

DKA法

代数方程式は重根も別々に考えると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が実軸に対称に分布することを防ぐためである(共役複素数になることを防ぐ).

ソースコード

上記のコードを含むVisual Studio 2010用のソースコードを以下に置く.

参考文献

  • 佐藤次男, 中村理一郎, ``よくわかる数値計算 アルゴリズムと誤差解析の実際'', 日刊工業新聞社, 2001.
  • 川上一郎, ``数値計算の基礎'', http://www7.ocn.ne.jp/~kawa1/

*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.

添付ファイル: filerootfinding.zip 742件 [詳細] filerootfinding.eq6.gif 444件 [詳細] filerootfinding.eq8.gif 453件 [詳細] filerootfinding.eq1.gif 464件 [詳細] filerootfinding.eq9.gif 449件 [詳細] filerootfinding.eq5.gif 403件 [詳細] filerootfinding.eq11.gif 420件 [詳細] filerootfinding.eq4.gif 433件 [詳細] filerootfinding.eq7.gif 414件 [詳細] filerootfinding.eq10.gif 421件 [詳細] filerootfinding.eq3.gif 468件 [詳細] filerootfinding.eq2.gif 416件 [詳細]

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