Newton Root Finder

Nie wymaga zakresu do poszukiwania miejsca zerowego, tylko punktu startowego. Wymaga za to obliczania wartości pochodnej, czyli także ciągłości funkcji. Metodę tą ciężko zastosować jeśli nie możemy podać funkcji na pochodną. Jest to także metoda iteracyjna. Krok iteracji ma postać:

$x_{1} = x_0 - \frac{f(x_0)}{f'(x_0)}$

Metoda ta prawie na pewno nie poradzi sobie z ekstremami funkcji - zacznie wokół nich oscylować. Może się zdarzyć, że zostanie znaleziony pierwiastek z poza zakresu, gdyż w jego kierunku iteracja zacznie podążać wzdłuż krzywej.

Poniższy kod stara się poszukuje miejsca zerowego dwa razy, obierając jako punkt startu a i b zakresu poszukiwań. Tym samym zwiększamy trochę szansę znalezienia pierwiastków. Ale ciągle nie radzi sobie z ekstremami.

Kod:

public class CentralDifferenceMethod
{
    private double m_h;
    private Func<double, double> m_func;

    public CentralDifferenceMethod()
        : base()
    {
    }
        
    public CentralDifferenceMethod(Func<double, double> a_func, 
        double a_h = Constants.DOUBLE_PRECISION * 2)
    {
        m_h = a_h;
        m_func = a_func;
    }

    public double Derive1(double a_x)
    {
        return (m_func(a_x + m_h) - m_func(a_x - m_h)) / (2 * m_h);
    }

    public double Derive2(double a_x)
    {
        return (m_func(a_x + m_h) - 2 * m_func(a_x) + 
            m_func(a_x - m_h)) / (m_h * m_h);
    }

    public static  Func<double, double> Derive1(
        Func<double, double> a_func, double a_h)
    {
        CentralDifferenceMethod cdm = 
            new CentralDifferenceMethod(a_func, a_h);

        return (x) =>
        {
            return cdm.Derive1(x);
        };
    }

    public static Func<double, double> Derive2(
        Func<double, double> a_func, double a_h)
    {
        CentralDifferenceMethod cdm = 
            new CentralDifferenceMethod(a_func, a_h);

        return (x) =>
        {
            return cdm.Derive2(x);
        };
    }

    public static Func<double, double> Derive2(
        Func<double, double> a_func, 
        Func<double, double> a_d1, double a_h)
    {
        CentralDifferenceMethod cdm = 
            new CentralDifferenceMethod(a_d1, a_h);

        return (x) =>
        {
            return cdm.Derive1(x);
        };
    }
}

public abstract class RootFinder
{
    public int MaxIterations;
    public double AbsoluteError;

    public RootFinder()
    {
        MaxIterations = Constants.ROOT_FINDER_MAXIMUM_ITERATIONS;
        AbsoluteError = Constants.ROOT_FINDER_ABSOLUTE_ERROR;
    }

    public static RootFinder Create(RootFindingMethod a_method)
    {
        switch (a_method)
        {
            case RootFindingMethod.Bisection: 
                return new BisectionRootFinder();
            case RootFindingMethod.ModifedRegulaFalsi: 
                return new ModifiedRegulaFalsiRootFinder();
            case RootFindingMethod.RegulaFalsi: 
                return new RegulaFalsiRootFinder();
            case RootFindingMethod.Secant: 
                return new SecantRootFinder();
            case RootFindingMethod.Newton: 
                return new NewtonRootFinder();
            case RootFindingMethod.Halley: 
                return new HalleyRootFinder();
            default: throw new NotImplementedException();
        }
    }

    public abstract double FindRoot(double a_a, double a_b, 
        Polynomial a_poly);
    public abstract double FindRoot(double a_a, double a_b, 
        Func<double, double> a_func);
    public abstract double FindRoot(double a_a, double a_b, 
        Func<double, double> a_func, Func<double, double> a_d1);
    public abstract double FindRoot(double a_a, double a_b, 
        Func<double, double> a_func, Func<double, double> a_d1, 
        Func<double, double> a_d2);
}

public abstract class DerivativeRootFinder : RootFinder
{
    protected struct RootState
    {
        public double Root;
        public bool Finded;

        public RootState(double a_x, bool a_finded = true)
        {
            Debug.Assert(a_x.IsNumber());

            Root = a_x;
            Finded = a_finded;
        }
    }

    protected double FindRootFromBothSide(double a_a, double a_b, 
        Func<double, double> a_func, Func<double, 
        RootState> a_root_finder)
    {
        Debug.Assert(a_a < a_b);
        Debug.Assert(a_a.IsNumber());
        Debug.Assert(a_b.IsNumber());

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

        if (Math.Abs(ya) < AbsoluteError)
            return a_a;
        if (Math.Abs(yb) < AbsoluteError)
            return a_b;

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

        RootState x1 = a_root_finder(a_a);
        if (x1.Finded && (x1.Root > a_a) && (x1.Root < a_b))
            return x1.Root;

        RootState x2 = a_root_finder(a_b);
        if (x2.Finded && (x2.Root > a_a) && (x2.Root < a_b))
            return x2.Root;

        return Double.PositiveInfinity;
    }
}

public class NewtonRootFinder : DerivativeRootFinder
{
    public override double FindRoot(double a_a, double a_b, 
        Func<double, double> a_func)
    {
        return FindRoot(a_a, a_b, a_func, 
            CentralDifferenceMethod.Derive1(a_func, AbsoluteError));
    }

    public override double FindRoot(double a_a, double a_b, 
        Func<double, double> a_func, Func<double, double> a_d1, 
        Func<double, double> a_d2)
    {
        return FindRoot(a_a, a_b, a_func, a_d1);
    }

    public override double FindRoot(double a_a, double a_b, 
        Func<double, double> a_func, Func<double, double> a_d1)
    {
        return FindRootFromBothSide(a_a, a_b, a_func, 
            (s) => FindRoot(s, a_func, a_d1));
    }

    private RootState FindRoot(double a_x, Func<double, double> a_func,
        Func<double, double> a_d1)
    {
        double prev_y = a_func(a_x);

        for (int i = 0; i < MaxIterations; i++)
        {
            double dy = a_d1(a_x);
            if (dy.IsAlmostEquals(0))
                return new RootState(a_x);

            a_x = a_x - prev_y / dy;
            double y = a_func(a_x);
            prev_y = y;

            if (Math.Abs(a_x) > AbsoluteError)
            {
                if (Math.Abs(y / a_x) < AbsoluteError)
                    return new RootState(a_x);
            }
            else if (Math.Abs(y) < AbsoluteError)
                return new RootState(a_x);
        }

        return new RootState(a_x, false);
    }

    public override double FindRoot(double a_a, double a_b, 
        Polynomial a_poly)
    {
        Polynomial d1 = a_poly.Differentiate();
        return FindRoot(a_a, a_b, (x) => a_poly.Evaluate(x),
            (x) => d1.Evaluate(x));
    }
}