Ostrowski's method

オストロフスキー法(Ostrowski's method)による求根

English edition is here.
Please use the ▶mirror site as well.

オストロフスキー法は根を求める方法の一つで、逆補間に多項式でなくパデ近似 (Padé Approximation) のように有理関数を用いている。他の方法に比べ安定性や収束速度の点で優れている。

逆補間に最も簡単な１次／１次の有理式 (1) を用いた場合には

… (1)

次式 (2) のように表せる。

… (2)

式 (2) の右辺を整理すると α δ - β γ の箇所が消えて次の関係が得られる。

… (3)

式 (3) に g(x) = 0 を代入して根の近似値が得られ、

根の近似値 x_n+3 は左辺の分子を t として次のようになる。

… (4)

… (5)

重根でない場合には根の近くでは 2次の収束を示す。

C プログラム例

list.1



/* ost - Solve the equation f(x) = 0 to within desired error value for x.
 *     Rev.3.14: Sep. 3, 2011
 *     Copyright(c) 1998-2011 Takayuki HOSODA, All rights reserved.
 *     http://www.finetune.jp/~lyuka/technote/fract/ost.html
 *
 * SYNOPSIS
 *     #include <math.h>
 *     #include <stdlib.h>
 *     #include <float.h>
 *
 *     double
 *     ost(double guess1, double guess2, double (*func)(), double epsilon, double *x);
 *
 * RETURN VALUE
 *     The ost() function returns the result code, and stores a root to the x
 *     when it found,
 *     0 : A root found.
 *     1 : An approximate solution within the specified epsilon found. |f(x)| <= epsilon
 *     2 : The search reached to the calculation limit. x_n+1 == x_n
 *     3 : The search concentrated in a local flat region of the function. f(x_n+2) == f(x_n+1)
 *     4 : The search reached to the calculation limit. |f(x_n+1)| > |f(x_n)|
 *     5 : No root found while specified iteration.
 *     6 : The search diverged.
 *     7 : Error. Bad guess or f(x) is constant.
 *     8 : Error. Bad argument, Iteration is negative.
 *     9 : Error. Bad argument, func() points NULL.
 *    10 : Error. Bad argument, x points NULL.
 */
 
#include <math.h>
#include <stdlib.h>
#include <float.h>
#ifdef DEBUG
#include <stdio.h>
#endif

#define MAXPVT             31
#define M2E31              2.147483648e9
#define EXP_THRESHOLD      1e-8
#define EXP_FACTOR         -1.5
#ifndef DBL_EPSILON
#define DBL_EPSILON        2.2204460492503131E-16
#endif

double nextc(a, b) double a, b; { return a + (b - a) * (1.0 + 14.0 * ((double)random() / M2E31)) / 16.0; }
double delta(a, b) double a, b; { return b == 0.0 ? fabs(a) : (fabs(a) - fabs(b)) / fabs(b); }
double minx(a, b, c, d, e, f) double a, b, c, d, e, f; {
    return fabs(f) < fabs(e) ? fabs(f) < fabs(d) ? c : a : fabs(e) < fabs(d) ? b : a;
}

int
ost(guess1, guess2, func, epsilon, iteration, x)
    double          guess1;
    double          guess2;
    double          (*func)(double);
    double          epsilon;
    int             iteration;
    double          *x;
{
    double          a, b, c, d, e, f, t;
    int             i;
    int             p = 0;
    double          q = 0.0;
    double          r0 = 0.0;
    static double   r1 = 0.0;
    static double   r2 = 0.0;
    int             xp =0;
    int             ft =0;

    if (func == NULL) return 10;
    if (x == NULL) return 9;
    if (0 >= iteration) { return 8; }
    a = guess1; b = guess2;
    b = (a == b) ? a == 0 ? 0.01 : a * 1.01 : b;
    if (0.0 == (d = (*func)(a))) { *x = a; return 0; }
    if (fabs(d) <= epsilon)  { *x = a; return 1; }
    if (0.0 == (e = (*func)(b))) { *x = b; return 0; }
    if (fabs(e) <= epsilon)  { *x = b; return 1; }
    c = (a + b) * 0.5; f = (*func)(c);
    r2 = c; r1 = 0.0;
    for (i = 0; i < iteration; i++) {
        if (0.0 == f) { *x = c; return 0; }
        if (fabs(f) <= epsilon)  { *x = c; return 1; }
        if (c == a) {
            if ((p > MAXPVT) || (ft >= 2)) { *x = minx(a, b, c, d, e, f); return 2;}
            ft++; p++; f = (*func)(c = nextc(a, b)); continue;
        }
        if (fabs(f) == INFINITY || f == NAN) {
            if (p > MAXPVT) { *x = minx(a, b, c, d, e, f); return 6; }
            i++; a = (a + b) * 0.5; d = (*func)(a);
            p++; c = nextc(a, b); f = (*func)(c); continue;
        }
        if (i >= 5 && (fabs(f) >= fabs(e) * 10.0)) {
            if (p > MAXPVT) { *x = minx(a, b, c, d, e, f); return 5; }
            if (fabs(d) > fabs(e) * 10.0) { i++; a = (a + b) * 0.5; d = (*func)(a); }
            p++; c = nextc(a, b); f = (*func)(c); continue;
        }
        if (ft >= 2 && xp != 3 && (fabs(f) > fabs(e)) && (EXP_THRESHOLD >= fabs(e)) && (i > p + 3)) {
            *x = minx(a, b, c, d, e, f); return 4;
        }
        if ((f == e) || (f == d)) {
            if (p > MAXPVT) {*x = c; return 7;}
            if (q == 0.0) q = -fabs(b - a);
            if (EXP_THRESHOLD < fabs(f)) {
                if (EXP_THRESHOLD > fabs(f + f - e - d)) {
                    b = (f * e < 0.0) ? nextc(b, c) : b + (q *= EXP_FACTOR); e = (*func)(b); p++; continue;
                }
                p++; c = nextc(a, b); f = (*func)(c); continue;
            } else {
                q = 0.0;
            }
            if (f == e && DBL_EPSILON  >= delta(c, a)) {*x = fabs(f) < fabs(d) ? c : a; return 2;}
            if (f == d && DBL_EPSILON  >= delta(c, b)) {*x = fabs(f) < fabs(e) ? c : b; return 2;}
            if (ft++ >= 2) {*x = minx(a, b, c, d, e, f); return 3; }
            p++; c = nextc(a, b); f = (*func)(c); continue;
        }
        t = e * (f - d) * (c - b) / ( d * (f - e) * (c - a));
        t  = (t * a - b) / (t - 1.0); /* Ostrowski's method */
        if (t == c) {
            if (ft >= 2) {*x = minx(a, b, c, d, e, f); return 2;}
            t += (ft == 1) ? -t * DBL_EPSILON : t * DBL_EPSILON; /* finetune */
        } else {
            if (fabs(d) > fabs(e) && fabs(e) > fabs(f)) ft = 0.0;
        }
        if ( 5.0 * fabs(fabs(e) - fabs(f)) < fabs(e) + fabs(f) && f * e < 0.0) {
            t = (f * b - e * c) / (f - e); /* Use regula falsi method */
        } else if (0 == xp) {
            r0 = r1; r1 -= (r1 / 0.53 + (c - b) / (b - a) - r2) * 0.5; r2 -= r0 * 0.5;  /* IIR filter */ 
            if ((0.5 < r2 && r2 < 1.0 && i >= 15 && EXP_THRESHOLD >= fabs(f)) && 0.01 > fabs(r1 / r2)) {
                xp = 3; t = c - (c - b) * (c - b) / (c - 2 * b + a); /* Use Aitken method instead */
            }
        } else {
            xp--;
        }
        a = b; b = c; c = t; d = e; e = f; f = (*func)(c);
    }
    *x = minx(a, b, c, d, e, f); return 5;
}

#ifdef DEBUG
#include <stdio.h>
#define ITERATION 100
#define EPSILON  2.2204460492503131E-16

double f0(x) double x; { return exp(x) - 3 * x * x; }
#ifdef DEBUG_MORE
double f1(x) double x; { return x * x * (x + 4) - 15.0; }
double f2(x) double x; { return sin(x) - x / 2.0; }
double f3(x) double x; { return exp(-x) + cos(x); }
double f4(x) double x; { return 10.0 * x * exp(-(x * x)) - 1.0 ; }
double f5(x) double x; { return atan(x) - x + 1.0; }
double f6(x) double x; { return tan(x) - 1.0 / x ; }
#endif
double f7(x) double x; { double t = x - 1.0; return (t * t * t * t * t) ; }

double
fs(x)
    double          x;
{
    double          t;
    double          p = 1e-5;

    t = x + 273.15;
    if (t < 0) return -p;
    return exp(((16.635794 + 2.433502 * log(t) + 0.00001673952 * t * t) - 0.02711193 * t)
               - 6096.9385 / t) * 100.0 - p;
}

int
main()
{
    double        root, a, b, eps;
    int           r;
    extern double f();

    printf("main:f0: exp(x) - 3 * x^2 == 0\n");
    r= ost(a=3.0, b=4.0, f0, eps=EPSILON, ITERATION, &root);
    printf("main:f0:[%.16g, %.16g], eps=%g,  f0(%.16g) = %.16g, r=%d\n", a, b, eps, root, f0(root), r);
    printf("main:f0: exp(x) - 3 * x * x == 0\n");
    r= ost(a=-1.0, b=0.0, f0, eps=EPSILON, ITERATION, &root);
    printf("main:f0:[%.16g, %.16g], eps=%g,  f0(%.16g) = %.16g, r=%d\n", a, b, eps, root, f0(root), r);
    printf("main:f0: exp(x) - 3 * x^2 == 0\n");
    r= ost(a=0.0, b=1.0, f0, eps=EPSILON, ITERATION, &root);
    printf("main:f0:[%.16g, %.16g], eps=%g,  f0(%.16g) = %.16g, r=%d\n", a, b, eps, root, f0(root), r);
    printf("main:fs: 'Water saturation pressure (Sonntag) == 1e-5'\n");
    r= ost(a=-273.15, b=100.0, fs, eps=EPSILON, ITERATION, &root);
    printf("main:fs:[%.16g, %.16g], eps=%g,  fs(%.16g) = %.16g, r=%d\n", a, b, eps, root, fs(root), r);
    printf("main:f7: (x - 1)^5 == 0\n");
    r= ost(a=0.0, b=1.5, f7, eps=1e-100, ITERATION, &root);
    printf("main:f7:[%.16g, %.16g], eps=%g,  f7(%.16g) = %.16g, r=%d\n", a, b, eps, root, f7(root), r);
    return 0;
}
#endif

list.1 を -DDEBUG オプション付きでコンパイルし実行した結果の例


main:f0: exp(x) - 3 * x^2 == 0
main:f0:[3, 4], eps=2.22045e-16,  f0(3.733079028632814) = 3.83026943495679e-15, r=2
main:f0: exp(x) - 3 * x * x == 0
main:f0:[-1, 0], eps=2.22045e-16,  f0(-0.4589622675369485) = -7.215365457891032e-17, r=1
main:f0: exp(x) - 3 * x^2 == 0
main:f0:[0, 1], eps=2.22045e-16,  f0(0.9100075724887091) = -2.244298497044994e-16, r=2
main:fs: 'Water saturation pressure (Sonntag) == 1e-5'
main:fs:[-273.15, 100], eps=2.22045e-16,  fs(-125.4179201811381) = -1.001192943654583e-18, r=1
main:f7: (x - 1)^5 == 0
main:f7:[0, 1.5], eps=1e-100,  f7(1) = 0, r=0

詳細は : Convergence details of the test functions

APPENDIX.A

Step response of the IIR filter used in the program shown in list.1

APPENDIX.B

オストロフスキーは根の近似に逆補間を用いる求根アルゴリズムの一つなので、ラグランジェ補完のような他の補完でも求根アルゴリズムとして使える。
もしそれの n=3 のを用いるとすれば、漸化式、
… (6)

が前述のオストロフスキーの漸化式 (4), (5)

の代わりに使える。
これは、list.1 中の一行
t = (t * a - b) / (t - 1.0);
を
t = a * e * f * (f - e) - b * d * f * (f - d) + c * d * e * (e - d); t /= (f - e) * (f - d) * (e - d);
で置き換えることによって試すことが出来る。

上記の方法では根の多重度が１のゼロの近傍ではほとんど２次（約1.8次) の収束を示すが、予測値が根の付近ではない場合には発散しがちであるため推奨出来ない。
私見を述べるなら、求根アルゴリズムとして最も大切なのは収束の次数ではなく、収束の安定さ、つまり滅多に発散せずに根を見つける能力である。

参考文献

Ostrowski Alexander M., 1966, Solutions of Equations and Systems of Equations, 2nd ed. (New york Academic Press), Chapter 12.
Ostrowski Alexander M., 1973, Solution of Equations in Euclidean and Bonach Spaces, Academic Press, New York, 3rd ed.