User Manual 4.17 Root-Finding Algorithms

Introduction

Scope

This section is about the root-finding algorithms for univariate real functions provided by the library. First we explain how to use the root-finding algorithms. Then a focus is made on the following root finding methods: Brent, Newton, Bisection and Müller.

Javadoc

All solvers are in the following package : [[[:Modèle:JavaDoc4.17]]/fr/cnes/sirius/patrius/math/analysis/solver/package-summary.html fr.cnes.sirius.patrius.math.analysis.solver]

Links

None as of now.

Useful Documents

None as of now.

Package Overview

The package fr.cnes.sirius.patrius.math.analysis.solvers contains all the solvers described here.

Features Description

About root finding

This library provides root-finding algorithms for real univariate functions. The function for which a root is searched has to be provided as a [[[:Modèle:JavaDoc4.17]]/fr/cnes/sirius/patrius/math/analysis/UnivariateFunction.html UnivariateFunction]. The root-finding algorithm is implemented through the [[[:Modèle:JavaDoc4.17]]/fr/cnes/sirius/patrius/math/analysis/solver/UnivariateSolver.html UnivariateSolver] interface.

A very generic example would be :

        UnivariateFunction f = <some function>;
        UnivariateSolver solver = <some solver>;
        int numberOfEvals = 100;    
        double startInterval = 5.;
        double endInterval = 10.;

        double root = solver.solve(numberOfEvals , f, startInterval , endInterval);

This means the solver with look for a root value of f in the interval [5. , 10.], and will have to do so in no more than 100 evaluation steps.

Test functions

We will demonstrate the use of root-finding algorithms with two test input functions : a periodic function and a linear function. These functions are defined as follows :
sinus function : SinFunction

public class SinFunction implements DifferentiableUnivariateFunction {

    public double value(double x) {
        return FastMath.sin(x);
    }

    public UnivariateFunction derivative() {
        return new UnivariateFunction() {
            public double value(double x) {
                return FastMath.cos(x);
            }
        };
    }

linear function : XMinus5Function

public class XMinus5Function implements DifferentiableUnivariateFunction {

    public double value(double x) {
        return x- 5;
    }

    public UnivariateFunction derivative() {
        return new UnivariateFunction() {
            public double value(double x) {
                return 1.0;
            }
        };
    }

Available solvers

The library provides around ten solvers, but the following have been more thoroughly tested :

• BrentSolver
• NewtonSolver
• BisectionSolver
• MullerSolver
• LaguerreSolver
• PolynomialRootsFinder

Brent method

TheBrent's method is a root-finding algorithm combining the bisection method, the secant method and inverse quadratic interpolation. It has the reliability of bisection but it can be as quick as some of the less reliable methods. The idea is to use the secant method or inverse quadratic interpolation if possible, because they converge faster, but to fall back to the more robust bisection method if necessary.

        UnivariateFunction f = new SinFunction ();
        double r;
        UnivariateSolver solver = new BrentSolver();

        r = solver.solve(100, f, 3, 4);

Result : r = PI

        UnivariateFunction f = new XMinus5Function ();
        double r;
        UnivariateSolver solver = new BrentSolver();

        r = solver.solve(100, f, 1, 2);

Result : NoBracketingException exception

behavior deteriorated in case of several roots in the interval :

        UnivariateFunction f = new SinFunction ();
        double r;
        UnivariateSolver solver = new BrentSolver();

        r = solver.solve(100, f, 1, 20);

Result : r = PI

        UnivariateFunction f = new SinFunction ();
        double r;
        UnivariateSolver solver = new BrentSolver();

        r = solver.solve(100, f, 2, 20);

Result : r = 3 PI

        UnivariateFunction f = new SinFunction ();
        double r;
        UnivariateSolver solver = new BrentSolver();

        r = solver.solve(100, f, 7, 20);

Result : NoBracketingException exception

Newton method

The Newton's method (also known as the Newton–Raphson method) is a method for finding successively better approximations to the roots of a real-valued function. The idea of the method is as follows: one starts with an initial guess which is reasonably close to the true root, then the function is approximated by its tangent line, and one computes the x-intercept of this tangent line. This x-intercept will typically be a better approximation to the function's root than the original guess, and the method can be iterated.

IMPORTANT : the Newton solver only works for differentiable functions, this is reflected in the interface and class inheritances.

        DifferentiableUnivariateFunction f = new XMinus5Function ();
        double r;
        DifferentiableUnivariateSolver solver = new NewtonSolver();

        r = solver.solve(100, f, 1, 6);

Result : r = 5

aberrant behavior :

        DifferentiableUnivariateFunction f = new SinFunction();
        double r;
        DifferentiableUnivariateSolver solver = new NewtonSolver();

        r = solver.solve(100, f, 1, 2);

Result : r =- 4 PI

        DifferentiableUnivariateFunction f = new SinFunction();
        double r;
        DifferentiableUnivariateSolver solver = new NewtonSolver();

        r = solver.solve(100, f, 1, 3);

Result : r = PI

        DifferentiableUnivariateFunction f = new XMinus5Function ();
        double r;
        DifferentiableUnivariateSolver solver = new NewtonSolver();

        r = solver.solve(100, f, 5.1, 20);

Result : r = 5

Bisection method

The bisection method is a root-finding method which repeatedly bisects an interval and then selects a subinterval in which a root must lie for further processing. It is a very simple and robust method, but it is also relatively slow.

        UnivariateFunction f = new XMinus5Function ();
        double r;
        UnivariateSolver solver = new BisectionSolver();

        r = solver.solve(100, f, 3, 4);

Result : r = PI

behavior deteriorated :

        UnivariateFunction f = new XMinus5Function ();
        double r;
        UnivariateSolver solver = new BisectionSolver();

        r = solver.solve(100, f, 5.1, 20);

Result : r = 20 (the upper bound of the interval)

Müller method

TheMüller's method is based on the secant method, which constructs at every iteration a line through two points on the graph of f. Instead, Müller's method uses three points, constructs the parabola through these three points, and takes the intersection of the x-axis with the parabola to be the next approximation.

        UnivariateFunction f = new XMinus5Function ();
        double r;
        UnivariateSolver solver = new MullerSolver();

        r = solver.solve(100, f, 1, 6);

Result : r = 5

Laguerre method

The Laguerre method is based on root finding of real coefficient polynomials. For reference, see “A First Course in Numerical Analysis, ISBN 048641454X, chapter 8.” In short, laguerre's method is global solver in the sense that it can start with any initial approximation and be able to solve all roots from that point. The algorithm requires a bracketing condition. When creating an instance of this solver, the relative accuracy, absolute accuracy and function value accuracy needs to be provided. With this instance, it is then possible to call the different “solveComplex” or “solveAllComplex” methods. The difference being that “solveComplex” methods find one complex root for the polynomial with the given coefficients, starting from the given initial value, while solveAllComplex methods return an array with all the complex roots of the polynomial. The polynomial is to be provided to the methods as a list of coefficients, in descending degree order. So, if for example we want to find the roots for 3x^2 – 2x + 4 = 0, the coefficients array should be [3, -2, 4]. It is also possible to configure the maximum number of iterations and the initial guess to start looking for solutions. Moreover, an array of initial guesses with their corresponding array of maximum iterations per guess can be given as well. If this is the case, the first initial point is checked. If no solution is found, the following initial points are checked. If after all the initial solution points the problem does not converge still, an error is raised. Possible errors:

• NullArgumentException if the coefficients of the polynomial are null
• NoDataException if the coefficients array is empty
• TooManyEvaluationsException if after the max number of iterations, the problem does not converge.

An example of code can be found below. In this example we try to find the complex solutions of a certain polynomial. Its coefficients are expressed in the “coefs” array. First, we check that using 0.0 as initial guess, the problem does not converge after 100 iterations. Then, we try to solve the problem again giving “0.0” and “-1.0” as initial guesses. The code should return all the complex roots of the polynomial:

        // Coefficients of the polynomial in descending degree order
        final double[] coefs = { -1277.3007688108385, -2784.4747201350438, -3484.5151185287687, -3025.129112243774, -1809.09448519026, -659.3527056579327, -192.0264504557372, -30.431131343092268, -4.143607944239528 };
        final LaguerreSolver laguerreSolver = new LaguerreSolver();
        
        // Check that with one initial point, no solution is obtained.
        try {
                final Complex[] rootsFail = laguerreSolver.solveAllComplex(coefs, 0.0, 100);
                // This part of the code should not be reached
        Assert.fail();
        } catch (TooManyEvaluationsException e) {
                // Check the message is the attended one.
                Assert.assertEquals(e.getMessage(), "illegal state: maximal count (100) exceeded: evaluations");
        }
        
        // Reconfigure the solver with multiple initial guesses:
        final Complex[] roots =
        laguerreSolver.solveAllComplex(coefs, new double[] { 0.0, -1.0 }, new int[] { 100, 100 });

        // Ref value of roots
        final Complex[] rootsRef = new Complex[] { new Complex(-1.0866082737090914, -0.5144910983744713),
        new Complex(-1.0866082737090914, 0.5144910983744713),
        new Complex(-0.19830816503367565, -1.1728862134932596),
        new Complex(-0.19830816503367568, 1.1728862134932596), new Complex(-1.3187864867933607, -3.034067652882213),
        new Complex(-1.3187864867933603, 3.034067652882214), new Complex(-1.0683543434675484, 3.553800385639665),
        new Complex(-1.068354343467548, -3.553800385639666) };
        
        // Verify value of roots
        for (int i = 0; i < roots.length; i++) {
                Assert.assertEquals(rootsRef[i].getReal(), roots[i].getReal(), Precision.DOUBLE_COMPARISON_EPSILON);
                Assert.assertEquals(rootsRef[i].getImaginary(), roots[i].getImaginary(), Precision.DOUBLE_COMPARISON_EPSILON);
        }

PolynomialRootsFinder method

The PolynomialRootsFinder method computes the roots of the polynomial p represented by an array of coefficients, starting with the coefficient of greater degree. A coefficient of 0 indicates that an intermediate power is not present in the equation. It can compute either just the real part of the roots, or the full complex solution.

In this solver, the roots of the polynomial are calculated by computing the eigenvalues of the companion matrix A:

        A = diag(ones(n-1, 1), -1);
        A(1,:) = -p(2:n+1)./p(1);
        r = eig(A);

The results produced are the exact eigenvalues of a matrix within roundoff error of the companion matrix, A. However, this does not mean that they are the exact roots of a polynomial whose coefficients are within roundoff error of those in p. The public method findRootsRealPolynomial() computes the real part of the roots of the polynomial, while findRootsComplexPolynomial() computes the full complex solution. To use either of these methods, it is only necessary to give as input the list of coefficients representing the polynomial (in decreasing degree order). So, if for example we want to find the roots for 3x^2 – 2x + 4 = 0, the coefficients array should be [3, -2, 4]. Possible errors:

• NullArgumentException if the coefficients of the polynomial are null
• NoDataException if the coefficients array is empty or the degree of the polynomial is lower than 1.

An example of code can be found below. In this example we try to find the complex solutions of a certain polynomial. Its coefficients are expressed in the “coefs” array.

        // Define coefficients
        final double[] coefficients = { 3, -2, 4 };
        
        // Compute roots
        double[] roots = PolynomialRootsFinder.findRootsRealPolynomial(coefficients);

Contents

Interfaces

The library defines the following interfaces related to solvers:

Interface	Summary	Javadoc
UnivariateSolver	Interface for a univariate solver	[[[:Modèle:JavaDoc4.17]]/fr/cnes/sirius/patrius/math/analysis/solver/UnivariateSolver.html ...]
UnivariateDifferentiableSolver	Interface for a differentiable univariate real solver	[[[:Modèle:JavaDoc4.17]]/fr/cnes/sirius/patrius/math/analysis/solver/UnivariateDifferentiableSolver.html ...]

Classes

This section is about the following classes related to solvers :

Class	Summary	Javadoc
BrentSolver	Brent solver implementation.	[[[:Modèle:JavaDoc4.17]]/fr/cnes/sirius/patrius/math/analysis/solver/BrentSolver.html ...]
NewtonRaphsonSolver	Newton-Raphson solver implementation.	[[[:Modèle:JavaDoc4.17]]/fr/cnes/sirius/patrius/math/analysis/solver/NewtonRaphsonSolver.html ...]
BisectionSolver	Bisection solver implementation.	[[[:Modèle:JavaDoc4.17]]/fr/cnes/sirius/patrius/math/analysis/solver/BisectionSolver.html ...]
MullerSolver	Müller solver implementation.	[[[:Modèle:JavaDoc4.17]]/fr/cnes/sirius/patrius/math/analysis/solver/MullerSolver.html ...]
LaguerreSolver	Laguerre solver implementation.	[[[:Modèle:JavaDoc4.17]]/fr/cnes/sirius/patrius/math/analysis/solver/LaguerreSolver.html ...]
PolynomialRootsFinder	PolynomialRootsFinder solver implementation.	[[[:Modèle:JavaDoc4.17]]/fr/cnes/sirius/patrius/math/analysis/solver/PolynomialRootsFinder.html ...]

Tips & Tricks

Summing up : the solvers appear to behave in unexpected ways when the initial conditions are not "ideal". Therefore It is advised, when using the solvers, to :

• avoid non-ideal initial conditions,
• read the solvers' javadoc carefully,
• always check if the solvers' results are appropriate.