Root Finder Modifed Regula Falsi

Metoda szybsza od dwóch pierwszej Bisection. Zawsze gwarantująca rozwiązanie w przeciwieństwie do metody Secant.

Wymaga przedziału startowego (a,b) takiego, że znaki f(a) i f(b) są różne, co oznacza, że w przedziale (a,b) istnieje miejsce zerowe. By tak się stało funkcja powinna być ciągła.

Jest to metoda rekurencyjna, w której stopniowo przybliżamy miejsce zerowe. W przedziale (a,b) mogą istnieć ekstrema nie będące miejsca miejscami zerowymi albo parzyste miejsca zerowe. Metoda ta nie zawsze gwarantuje nam znalezienie rozwiązania.

W dowolnym kroku iteracji mamy przedział $(x_0, x_1)$, podobnie jak w metodzie Secant wyznaczamy linię przebiegającą przez punkty $(x_0,f(x_0))$ i $(x_1,f(x_1))$. Znajdujemy jej miejsce zerowe:

$x_2=x_1-f(x_1)\frac{x_1-x_0}{f(x_1)-f(x_0)}=\frac{f(y_1)x_0-f(y_0)x_1}{f(x_1)-f(x_0)}$

Teraz wybieramy nowy zakres tak by znaki funkcji na końcach zakresu miały przeciwny znak. Cały proces powtarzamy.

Kod:

public override double FindRoot(
    double a_a, double a_b, Func a_func)
{
    Debug.Assert(a_a < a_b);

    double ya = a_func(a_a);
    double yb = a_func(a_b);

    if (Math.Abs(ya) < Constants.ROOT_FINDER_PRECISION)
        return a_a;

    if (Math.Abs(yb) < Constants.ROOT_FINDER_PRECISION)
        return a_b;

    if (Math.Sign(ya) * Math.Sign(yb) > 0.0)
        return Double.PositiveInfinity;

    double x = 0;

    for (int i = 0; i < Constants.ROOT_FINDER_MAXIMUM_ITERATIONS; i++)
    {
        x = (yb * a_a - ya * a_b) / (yb - ya);
        double y = a_func(x);

        if (Math.Abs(x) > Constants.ROOT_FINDER_PRECISION)
        {
            if (Math.Abs(y / x) < Constants.ROOT_FINDER_PRECISION)
                return x;
        }
        else if (Math.Abs(y) < Constants.ROOT_FINDER_PRECISION)
            return x;

        if ((Math.Sign(ya) * Math.Sign(y)) < 0)
        {
            a_b = x;
            yb = y;
        }
        else
        {
            a_a = x;
            ya = y;
        }
    }

    return x;
}

Metoda ta nie zawsze gwarantuje nam zbieżność. Istnieje jej modyfikacja która zawsze gwarantuje nam zbieżność. To dzielenie przez 2 ma uzasadnienie teorii.

public override double FindRoot(
    double a_a, double a_b, Func<double, double> a_func)
{
    Debug.Assert(a_a < a_b);

    double ya = a_func(a_a);
    double yb = a_func(a_b);

    if (Math.Sign(ya) * Math.Sign(yb) > 0.0)
        return Double.PositiveInfinity;

    if (Math.Abs(ya) < Constants.ROOT_FINDER_PRECISION)
        return a_a;

    if (Math.Abs(yb) < Constants.ROOT_FINDER_PRECISION)
        return a_b;

    double prev_y = ya;
    double x = 0;

    for (int i = 0; i < Constants.ROOT_FINDER_MAXIMUM_ITERATIONS; i++)
    {
        x = (yb * a_a - ya * a_b) / (yb - ya);
        double y = a_func(x);

        if (Math.Abs(x) > Constants.ROOT_FINDER_PRECISION)
        {
            if (Math.Abs(y / x) < Constants.ROOT_FINDER_PRECISION)
                return x;
        }
        else if (Math.Abs(y) < Constants.ROOT_FINDER_PRECISION)
            return x;

        if ((Math.Sign(ya) * Math.Sign(y)) < 0)
        {
            a_b = x;
            yb = y;
            if ((prev_y * y) > 0)
                ya /= 2;
        }
        else
        {
            a_a = x;
            ya = y;
            if ((prev_y * y) > 0)
                yb /= 2;
        }

        prev_y = y;
    }

    return x;
}