Class RootFinder<P>

java.lang.Object
org.bzdev.math.RootFinder<P>
Direct Known Subclasses:
RootFinder.Brent, RootFinder.Halley, RootFinder.Newton

public abstract class RootFinder<P> extends Object
Root finder.

Solves the equation y = f(x,p) for x, where x and y have a type of double and p has a type set by a type parameter and represents a set of parameters that do not vary. There are Three implementations: RootFinder.Brent, RootFinder.Newton and RootFinder.Halley that implement Brent's, Newton's and Halley's method respectively. If the type parameter is not used, one should not be provided. Otherwise it is a user-defined type used in computing the values of f and its derivatives with respect to x.

Once a RootFinder rf is constructed, calling rf.findRoot(guess) finds a root by solving 0 = f(x, p). To solve the equation y = f(x, p), call solve(y, guess) for the Newton and Halley case, or solve(y, lower, upper) for the Brent case, and the value of x will be returned. The argument "guess" is an initial guess and the arguments "lower" and "upper" are values of x bracketing the value giving the solution, where f(lower, p) - y and f(upper, p) - y have opposite signs.

Parameters allow a single root finder to use a number of related functions without the need to allocate a new root finder for each case. If parameters are used, one will define a class to hold the parameters. For example,


    class Parameters {
        double p1;
        double p2;
        public Parameters(double x1, double y1) {
           p1 = x1;
           p2 = x2;
        }
    }
 
Then create the RootFinder:

    Parameters parameters = new Parameters(1.0, 2.0);
    RootFinder<Parameters> rf = new RootFinder.Newton<Parameters>(parameters)
       {
          public double function(double x) {
             Parameters p = getParameters();
             double p1 = p.x1;
             double p2 = p.x2;
             double result;
             ...
             return result;
          }
          public double firstDerivative(double x) {
             Parameters p = getParameters();
             double p1 = p.x1;
             double p2 = p.x2;
             double result;
             ...
             return result;
          }
       };
 
To change the parameters, either create a new instance parameters2 and call rf.setParameters(parameters2) or call rf.getParameters() and manipulate the values (e.g., by setting parameters.p1 and parameters.p2).

The criteria for convergence can be changed by calling setEpsilon. Limits on the number of iteration steps and some control over those can be set by calling setLimits.

  • Nested Class Summary

    Nested Classes
    Modifier and Type
    Class
    Description
    static class 
    RootFinder class using Brent's method.
    static class 
    RootFinder exception class.
    static class 
    RootFinder class using Halley's method.
    static class 
    RootFinder class using Newton's method.
  • Field Summary

    Fields
    Modifier and Type
    Field
    Description
    protected int
    The iteration limit.
  • Constructor Summary

    Constructors
    Modifier
    Constructor
    Description
    protected
    Constructor.
    protected
    RootFinder(P parameters)
    Constructor given parameters.
  • Method Summary

    Modifier and Type
    Method
    Description
    abstract double
    ferror(double x)
    The error in the value returned by function(x).
    final double
    findRoot(double... initialArgs)
    Find a root.
    abstract double
    function(double x)
    The function f in the equation y = f(x, p).
    double
    Get the current value of epsilon The value of epsilon determines when the algorithm terminates normally.
    Get a RootFinder's parameters.
    static double
    refineSolution(DoubleUnaryOperator f, DoubleUnaryOperator df, double y, double guess)
    Use Newton's method to refine a solution to the equation f(x) - y = 0.
    void
    setEpsilon(double epsilon)
    Set the current value of epsilon.
    void
    setEpsilon(double epsilon, boolean relative)
    Set the current value of epsilon, indicating if it is a relative or absolute value.
    void
    setLimit(int maxLimit)
    Set the iterations limit.
    void
    setParameters(P parameters)
    Set a RootFinder's parameters.
    abstract double
    solve(double y, double... initialArgs)
    Solve an equation.
    static int
    solveBezier(double[] eqn, int offset, int n, double[] res)
    Find the zeros of a sum of Bernstein polynomial of degree n.
    static int
    solveCubic(double[] eqn)
    Solves the cubic equation ax3 + bx2 + cx + d and places the real roots into its argument array, which is used for both input and output.
    static int
    solveCubic(double[] eqn, double[] res)
    Solves the cubic equation ax3 + bx2 + cx + d = 0.
    static int
    solvePolynomial(double[] eqn, int n, double[] res)
    Find the zeros of a polynomial of degree n.
    static int
    Find the zeros of a BezierPolynomial.
    static int
    solvePolynomial(Polynomial p, double[] res)
    Find the zeros of a Polynomial.
    static int
    solveQuadratic(double[] eqn)
    Solves the quadratic equation ax2 + bx + c and places the real roots into its argument array, which is used for both input and output.
    static int
    solveQuadratic(double[] eqn, double[] res)
    Solves the quadratic equation ax2 + bx + c.
    static int
    solveQuartic(double[] eqn)
    Solve the quartic equation ax4 + bx3 + cx2 + dx + e = 0 and store the results in the argument array.

    Methods inherited from class java.lang.Object

    clone, equals, finalize, getClass, hashCode, notify, notifyAll, toString, wait, wait, wait
  • Field Details

    • maxLimit

      protected int maxLimit
      The iteration limit.
  • Constructor Details

    • RootFinder

      protected RootFinder()
      Constructor.
    • RootFinder

      protected RootFinder(P parameters)
      Constructor given parameters.
      Parameters:
      parameters - the parameters
  • Method Details

    • setParameters

      public void setParameters(P parameters)
      Set a RootFinder's parameters. Parameters are used to provide additional values for computing. root-finder's function and their derivatives.
      Parameters:
      parameters - the parameters
    • getParameters

      public P getParameters()
      Get a RootFinder's parameters. Parameters are used to provide additional values for computing. root-finder's function and their derivatives.
      Returns:
      an instance of the class representing a root finder's parameters (this will be the same instance passed to setParameters)
    • findRoot

      public final double findRoot(double... initialArgs) throws RootFinder.ConvergenceException
      Find a root. Starting from an initial guess, findRoot returns the value of x that satisfies f(x, p) = 0, where p represents the current parameters. The implementation simply calls solve with the parameter y set to 0, and with the initial arguments that solve requires (which depend on the root-finding algorithm being used).
      Parameters:
      initialArgs - the initial arguments required by a subclass
      Returns:
      the root
      Throws:
      RootFinder.ConvergenceException - the method failed to converge
      MathException - an error occurred calling the function f or one of its derivatives.
    • solve

      public abstract double solve(double y, double... initialArgs) throws RootFinder.ConvergenceException
      Solve an equation. Starting from an initial guess, findRoot returns the value of x that satisfies f(x, p) = y, where p represents the current parameters.
      Parameters:
      y - the desired value of f(x, p)
      initialArgs - the initial arguments required by a subclass
      Returns:
      the value of x that satisfies f(x, p) = y
      Throws:
      RootFinder.ConvergenceException - the method failed
      MathException - an error occurred calling the function f or one of its derivatives.
    • setLimit

      public void setLimit(int maxLimit)
      Set the iterations limit. The limit is a bound on the maximum number of iterations.
      Parameters:
      maxLimit - the maximum-iteration limit
    • getEpsilon

      public double getEpsilon()
      Get the current value of epsilon The value of epsilon determines when the algorithm terminates normally. If the current value of x is in the range [-1.0, 1.0] the requirement for convergence is that the change in the value of x from the previous iteration is less than epsilon and that the corresponding value of y is within epsilon of its desired value. Otherwise the requirement for convergence is that the change in the value of x, with that change divided by x, is less than epsilon and that the corresponding value of y is within epsilon of its desired value.
      Returns:
      the value of epsilon
    • setEpsilon

      public void setEpsilon(double epsilon) throws IllegalArgumentException
      Set the current value of epsilon. The value of epsilon determines when the algorithm terminates normally. The error associated with evaluating the function will be set to epsilon.
      Parameters:
      epsilon - the value of epsilon
      Throws:
      IllegalArgumentException - the argument was not a positive real number
    • setEpsilon

      public void setEpsilon(double epsilon, boolean relative) throws IllegalArgumentException
      Set the current value of epsilon, indicating if it is a relative or absolute value. When absolute, the error in evaluating the function will be epsilon when the default implementation of ferror(double) is used. When relative, the error will be the value of the function after that value is multiplied by epsilon.
      Parameters:
      epsilon - the value of epsilon
      relative - true if the value is relative; false if it is absolute.
      Throws:
      IllegalArgumentException - the first argument was not a positive real number
    • function

      public abstract double function(double x)
      The function f in the equation y = f(x, p). The function f(x, p) is used with x varied and with the parameters p constant.
      Parameters:
      x - the argument of the function
      Returns:
      the value of the function f
      See Also:
    • ferror

      public abstract double ferror(double x)
      The error in the value returned by function(x). The default implementation returns the same value as getEpsilon(). Callers are encouraged to override this method.
      Parameters:
      x - the argument passed to function(double)
      Returns:
      the error for the value function(x)
    • solveQuadratic

      public static int solveQuadratic(double[] eqn)
      Solves the quadratic equation ax2 + bx + c and places the real roots into its argument array, which is used for both input and output.

      This method duplicates the behavior of the method QuadCurve2D.solveQuadratic(double[]). It is provided to break a dependency on the jdk.desktop module.

      Parameters:
      eqn - the input/output array containing the roots (starting at index 0) for output and the coefficients ordered so that eqn[0] = c, eqn[1] = b, and eqn[2] = a.
      Returns:
      the number of roots; -1 if the coefficients a and b are zero
    • solveQuadratic

      public static int solveQuadratic(double[] eqn, double[] res)
      Solves the quadratic equation ax2 + bx + c. The length of the array res should be large enough to hold all of the roots (a length of at least 2 should be used).

      If multiple roots have the same value, only one is shown. If roots are sufficiently close that they cannot be distinguished due to floating-point limitations, only one root may be shown. This method duplicates the calling convention used by the method QuadCurve2D.solveQuadratic(double[],double[]). It was provided to break a dependency on the jdk.desktop module.

      Parameters:
      eqn - the input/output array containing the roots (starting at index 0) for output and the coefficients ordered so that eqn[0] = c, eqn[1] = b, and eqn[2] = a.
      res - the output array containing the roots (starting at index 0)
      Returns:
      the number of roots; -1 if the coefficients a and b are zero
    • solveCubic

      public static int solveCubic(double[] eqn)
      Solves the cubic equation ax3 + bx2 + cx + d and places the real roots into its argument array, which is used for both input and output.

      If multiple roots have the same value, only one of those roots is shown. If roots are sufficiently close that they cannot be distinguished due to floating-point limitations, only one root may be shown even though multiple distinct roots actually exist. This method duplicates the calling convention used by the method CubicCurve2D.solveCubic(double[]). It was provided to break a dependency on the jdk.desktop module, but experience with it indicates that is is more accurate than the implementation of CubicCurve2D.solveCubic(double[]) provided in Java 11.0.1.

      A cubic function f(x) = ax3 + bx2 + cx + d has an inflection point at -b/(3a). The number of critical points is dependent on the value Δ0 = b2 - 3ac:

      • if Δ is negative, there are no critical points and the cubic function is strictly monotonic. The implementation uses Newton's method, with the inflection point as an initial guess.
      • if Δ is zero, there is one critical point, which is also an inflection point. There is only one real root, which is a triple root if the inflection point is a root, or a single root if the inflection point is not a root. In the later case, the root's value is (-b - (27ad - b3)1/3 / (3a).
      • if Δ is positive, there are two critical points with the inflection point between these two. The critical points are $(-b-\sqrt{\Delta_0})/(3a)$ and $(-b + sqrt{\Delta_0})/(3a)$. If one of these is a root r (both cannot be roots), it is a double root, and the other root has the value -(2ar + b)/a. Otherwise, let x1 be the lower of these two critical points and let x2 be the higher of these two values.
        • If f(x1) and f″(x1) have opposite signs, there is a real root below x1. One can find the root using Newton's method with an initial guess $x_g = x_1 - (\frac32)^m\sqrt{-2f(x_1)/f\prime(x_1)}$, , where m is the smallest non-negative integer such that f(xg) and f(x1 have opposite signs.
        • If f(x1) and f(x1) have opposite signs, there is a root between these two values. One can find the root using Newton's method with a lower and upper bound and with the inflection point as an initial guess. When upper and lower bounds are provided, the RootFinder.Newton class will temporarily switch to Brent's method when Newton's method's next value is out of range, and will tighten the upper and lower bounds as new values are tried.
        • If f(x2) and f″(x2) have opposite signs, there is a real root above x2. One can find the root using Newton's method with an initial guess $x_g = x_2 - (\frac32)^m\sqrt{-2f(x_2)/f\prime(x_2)}$, where m is the smallest non-negative integer such that f(xg) and f(x2 have opposite signs.
      Parameters:
      eqn - the input/output array containing the roots (starting at index 0) for output and the coefficients ordered so that eqn[0] = d, eqn[1] = c, eqn[2] = b, and eqn[3] = a.
      Returns:
      the number of roots; -1 if the coefficients a, b, and c are zero
    • refineSolution

      public static double refineSolution(DoubleUnaryOperator f, DoubleUnaryOperator df, double y, double guess)
      Use Newton's method to refine a solution to the equation f(x) - y = 0. Newton's method uses a guess that is improved by approximating f(x) by x0 + f'(x0)(x - x0) and solving the linear equation y = x0 + f'(x0)(x - x0), iterating until a convergence criteria is satisfied or until the implementation detects that it will not converge. This implementation assumes the convergence criteria has been satisfied and iterates as long as the solution can be improved.

      The algorithm starts by setting x1 to the value of a guess that is provided by one of this methods arguments.

      A sequence of values {x1, ... xn, ...} is generated as follows: For xn, if f(xn) - y = 0, then xn+1 = x. Otherwise if f(xn) - y != 0, Newton's method is used to determine a value w = xn = (f(xn)-y)/f'(x1)) unless f'(x1) = 0, in which case w = xn. If the sign if f(xn) - y and f(w) - y differ, a point z = x1 + (u - x1) (f(x1 - y) / (f(x1) - f(w)) is also considered. Note that z is an interpolation point between x1 and w. The value of |f(x) - y| is compared at x1, w, and (if calculated) z and the value of xn+1 is set to the one which provides the minimal value of |f(x)-y|.

      This sequence converges when xn+1 = xn, and that must occur because there are a finite number of double-precision numbers between f(x1) and 0.

      The final step in the algorithm is to compute f(x) at adjacent floating-point numbers and pick the one that provides the minimum value for |f(x) - y|. The initial guess will be used in this final step if the argument df is null.

      Parameters:
      f - the function
      df - the function's derivative; null if only the final step is used
      y - the value for f(x)
      guess - an initial guess (assumed to be very close to the solution)
      Returns:
      a value x such that f(x) = y to within floating-point errors
    • solveCubic

      public static int solveCubic(double[] eqn, double[] res)
      Solves the cubic equation ax3 + bx2 + cx + d = 0. The length of the array res should be large enough to hold all of the real roots (a length of at least 3 should be used unless one knows a priori that there are less than 3 real roots).

      If multiple real roots have the same value, only one of those roots is shown. If roots are sufficiently close that they cannot be distinguished due to floating-point limitations, only one root may be shown even though multiple distinct roots actually exist. This method duplicates the calling convention used by the method CubicCurve2D.solveCubic(double[],double[]). It was provided to break a dependency on the jdk.desktop module, but experience with it indicates that is is more accurate than the implementation of CubicCurve2D.solveCubic(double[],double[]) provided in Java 11.0.1.

      A cubic function f(x) = ax3 + bx2 + cx + d has an inflection point at -b/(3a). The number of critical points is dependent on the value Δ0 = b2 - 3ac:

      • if Δ is negative, there are no critical points and the cubic function is strictly monotonic. The implementation uses Newton's method, with the inflection point as an initial guess.
      • if Δ is zero, there is one critical point, which is also an inflection point. There is only one real root, which is a triple root if the inflection point is a root, or a single root if the inflection point is not a root. In the later case, the root's value is (-b - (27ad - b3)1/3 / (3a).
      • if Δ is positive, there are two critical points with the inflection point between these two. The critical points are $(-b - \sqrt{\Delta_0})/(3a)$ and $(-b + \sqrt{\Delta_0})/(3a)$. If one of these is a root r (both cannot be roots), it is a double root, and the other root has the value -(2ar + b)/a. Otherwise, let x1 be the lower of these two critical points and let x2 be the higher of these two values.
        • If f(x1) and f″(x1) have opposite signs, there is a real root below x1. One can find the root using Newton's method with an initial guess $x_g = x_1 - (\frac32)^m\sqrt{-2f(x_1)/f\prime(x_1)}$, where m is the smallest non-negative integer such that f(xg) and f(x1 have opposite signs.
        • If f(x1) and f(x1) have opposite signs, there is a root between these two values. One can find the root using Newton's method with a lower and upper bound and with the inflection point as an initial guess. When upper and lower bounds are provided, the RootFinder.Newton class will temporarily switch to Brent's method when Newton's method's next value is out of range, and will tighten the upper and lower bounds as new values are tried.
        • If f(x2) and f″(x2) have opposite signs, there is a real root above x2. One can find the root using Newton's method with an initial guess $x_g = x_2 - (\frac32)^m\sqrt{-2f(x_2)/f\prime(x_2)}$, xg = x2 - (3/2)msqrt(-2f″(x2)/f(x2), where m is the smallest non-negative integer such that f(xg) and f(x2 have opposite signs.
      Parameters:
      eqn - the input/output array containing the roots (starting at index 0) for output and the coefficients ordered so that eqn[0] = d, eqn[1] = c, eqn[2] = b, and eqn[3] = a.
      res - the output array containing the roots (starting at index 0) in numerical order
      Returns:
      the number of roots; -1 if the coefficients a, b, and c are zero
    • solvePolynomial

      public static int solvePolynomial(Polynomial p, double[] res)
      Find the zeros of a Polynomial.
      Parameters:
      p - the polynomial
      res - an array to hold the results
      Returns:
      the number of roots; -1 if the degree of the polynomial is zero.
      Throws:
      IllegalArgumentException - if the array argument is too short
    • solvePolynomial

      public static int solvePolynomial(double[] eqn, int n, double[] res)
      Find the zeros of a polynomial of degree n.
      Parameters:
      eqn - the coefficients of the polynomial for input, ordered so that eqn[i] provides the coefficient for xi where x is the polynomial's indeterminate
      n - the degree of the polynomial
      res - the roots that were found
      Returns:
      the number of roots; -1 if the degree of the polynomial is zero.
      Throws:
      IllegalArgumentException - n is negative or the array arguments are too short
    • solveQuartic

      public static int solveQuartic(double[] eqn)
      Solve the quartic equation ax4 + bx3 + cx2 + dx + e = 0 and store the results in the argument array. The algorithm is based on the description in the Wikipedia article on the quartic function.

      A number of special cases are handled explicitly and for coefficients that are integer values, a computation of the discriminant is performed using two's complement arithmetic with long integers. This is equivalent to modular arithmetic (modulo 264) and allows one to reliably determine if the discriminant is zero, and to improve the accuracy of determining if the discriminant is positive, negative, or zero. Also, the real-valued form of the discriminant uses Kahan's additional algorithm to improve the accuracy of the sum.The same procedure is used for some other quantities as well. These are combined to determine the number of roots based on whether these quantities are positive, zero, or negative. Cases that are not covered should not occur unless as the result of floating point errors, in which case a numerical method is used instead. To further reduce floating-point errors, Newton's method and some floating-point adjustments are made to improve the accuracy of the solutions. In addition, for cases where there are known to be multiple roots, the critical points are computed by solving a cubic equation and those solutions are used when they are also solutions. If floating-point errors split a single multiple root of the polynomial into multiple zeros, the critical points are used to eliminate these spurious solutions.

      Parameters:
      eqn - an array containing the coefficients of the polynomial, where eqn[0] = e, eqn[1] = d, eqn[2] = c, eqn[3] = d, and eqn[4] = e (the array length must be at least 5), also used as the output array
      Returns:
      the number of real-valued solutions; -1 if a, b, c, and d are 0
    • solvePolynomial

      public static int solvePolynomial(BezierPolynomial p, double[] res)
      Find the zeros of a BezierPolynomial.
      Parameters:
      p - the polynomial
      res - an array to hold the results
      Returns:
      the number of roots; -1 if the degree of the polynomial is zero.
      Throws:
      IllegalArgumentException - if the array argument is too short
    • solveBezier

      public static int solveBezier(double[] eqn, int offset, int n, double[] res) throws IllegalArgumentException
      Find the zeros of a sum of Bernstein polynomial of degree n. The name of the method, solveBezier, reflects the use of Bernstein polynomials for Bézier curves. The roots returned are in the range [0.0, 1.0], the range appropriate for Bézier curves.
      Parameters:
      eqn - the coefficients of the polynomial for input, where eqn[i] is the coefficient for the Bernstein polynomial Bi,n(t)
      offset - the offset into the array eqn at which the coefficients start
      n - the degree of the polynomial
      res - the roots that were found
      Returns:
      the number of roots; -1 if the degree of the polynomial is zero.
      Throws:
      IllegalArgumentException - n is negative or the array arguments are too short