Secant Root Finder

Kolejna metoda na znajdowanie pierwiastków równania. Do rozpoczęcia wymagane są dwa punkty a i b, które powinny być dostatecznie blisko miejsca zerowego. Nie jest wymagane by miejsce zerowe znajdowało się wewnątrz (a,b), ale wtedy wzrasta ryzyko, że metoda może nie znaleźć miejsca zerowego. Podobnie dzieje się jeśli w zakresie (a,b) istnieje ekstremum nie będące miejscem zerowym. Tak więc możemy powiedzieć, że zakres (a,b) jest mniejszy tym lepiej.

Przybliżanie następuje rekurencyjnie. W dowolnym momencie iteracji mamy dwa punkty $y_0=f(x_0)$ i $y_1=f(y_1)$. Wyznaczamy równanie linii przez nie przechodzące:

$y=\displaystyle\frac{f(x_1)-f(x_0)}{x_1-x_0}(x-x_1) + f(x_1)$

Wyznaczym punkt przecięcia się naszej funkcji z linią:

$0=\displaystyle\frac{f(x_1)-f(x_0)}{x_1-x_0}(x-x_1) + f(x_1)$

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

Bierzemy nowy przedział $(x_1, x_2)$ i powtarzamy iterację. Ponieważ metoda nie zawsze znajduje rozwiązanie powinniśmy wyznaczyć górną ilość iteracji. Normalnie metoda powinna zakończyć działanie po znalezieniu takie x dla, którego f(x) jest w granicach precyzji.

Kod:

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.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 = a_a - ya * (a_a - a_b) / (ya - yb);

        Debug.Assert(x.IsNumber());

        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;

        yb = ya;
        a_b = a_a;

        ya = y;
        a_a = x;

        if (ya == yb)
            return Double.PositiveInfinity;
    }

    return x;
}