非線形方程式の根

  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;
}

  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;
}

  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;
}

  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;
}

  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;
}