BrentSolver.java

  1. /*
  2.  * Licensed to the Apache Software Foundation (ASF) under one or more
  3.  * contributor license agreements.  See the NOTICE file distributed with
  4.  * this work for additional information regarding copyright ownership.
  5.  * The ASF licenses this file to You under the Apache License, Version 2.0
  6.  * (the "License"); you may not use this file except in compliance with
  7.  * the License.  You may obtain a copy of the License at
  8.  *
  9.  *      http://www.apache.org/licenses/LICENSE-2.0
  10.  *
  11.  * Unless required by applicable law or agreed to in writing, software
  12.  * distributed under the License is distributed on an "AS IS" BASIS,
  13.  * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
  14.  * See the License for the specific language governing permissions and
  15.  * limitations under the License.
  16.  */
  17. package org.apache.commons.numbers.rootfinder;

  19. import java.util.function.DoubleUnaryOperator;

  21. /**
  22.  * This class implements the <a href="http://mathworld.wolfram.com/BrentsMethod.html">
  23.  * Brent algorithm</a> for finding zeros of real univariate functions.
  24.  * The function should be continuous but not necessarily smooth.
  25.  * The {@code solve} method returns a zero {@code x} of the function {@code f}
  26.  * in the given interval {@code [a, b]} to within a tolerance
  27.  * {@code 2 eps abs(x) + t} where {@code eps} is the relative accuracy and
  28.  * {@code t} is the absolute accuracy.
  29.  * <p>The given interval must bracket the root.</p>
  30.  * <p>
  31.  *  The reference implementation is given in chapter 4 of
  32.  *  <blockquote>
  33.  *   <b>Algorithms for Minimization Without Derivatives</b>,
  34.  *   <em>Richard P. Brent</em>,
  35.  *   Dover, 2002
  36.  *  </blockquote>
  37.  */
  38. public class BrentSolver {
  39.     /** Relative accuracy. */
  40.     private final double relativeAccuracy;
  41.     /** Absolute accuracy. */
  42.     private final double absoluteAccuracy;
  43.     /** Function accuracy. */
  44.     private final double functionValueAccuracy;

  46.     /**
  47.      * Construct a solver.
  48.      *
  49.      * @param relativeAccuracy Relative accuracy.
  50.      * @param absoluteAccuracy Absolute accuracy.
  51.      * @param functionValueAccuracy Function value accuracy.
  52.      */
  53.     public BrentSolver(double relativeAccuracy,
  54.                        double absoluteAccuracy,
  55.                        double functionValueAccuracy) {
  56.         this.relativeAccuracy = relativeAccuracy;
  57.         this.absoluteAccuracy = absoluteAccuracy;
  58.         this.functionValueAccuracy = functionValueAccuracy;
  59.     }

  61.     /**
  62.      * Search the function's zero within the given interval.
  63.      *
  64.      * @param func Function to solve.
  65.      * @param min Lower bound.
  66.      * @param max Upper bound.
  67.      * @return the root.
  68.      * @throws IllegalArgumentException if {@code min > max}.
  69.      * @throws IllegalArgumentException if the given interval does
  70.      * not bracket the root.
  71.      */
  72.     public double findRoot(DoubleUnaryOperator func,
  73.                            double min,
  74.                            double max) {
  75.         // Avoid overflow computing the initial value: 0.5 * (min + max)
  76.         // Note: This sum is invalid if min == max == Double.MIN_VALUE
  77.         // so detect this edge case. It will raise a bracketing exception
  78.         // if min is not the root within the configured function accuracy;
  79.         // otherwise min is returned.
  80.         final double initial = min == max ? min : 0.5 * min + 0.5 * max;
  81.         return findRoot(func, min, initial, max);
  82.     }

  84.     /**
  85.      * Search the function's zero within the given interval,
  86.      * starting from the given estimate.
  87.      *
  88.      * @param func Function to solve.
  89.      * @param min Lower bound.
  90.      * @param initial Initial guess.
  91.      * @param max Upper bound.
  92.      * @return the root.
  93.      * @throws IllegalArgumentException if {@code min > max} or
  94.      * {@code initial} is not in the {@code [min, max]} interval.
  95.      * @throws IllegalArgumentException if the given interval does
  96.      * not bracket the root.
  97.      */
  98.     public double findRoot(DoubleUnaryOperator func,
  99.                            double min,
  100.                            double initial,
  101.                            double max) {
  102.         if (min > max) {
  103.             throw new SolverException(SolverException.TOO_LARGE, min, max);
  104.         }
  105.         if (initial < min ||
  106.             initial > max) {
  107.             throw new SolverException(SolverException.OUT_OF_RANGE, initial, min, max);
  108.         }

  110.         // Return the initial guess if it is good enough.
  111.         final double yInitial = func.applyAsDouble(initial);
  112.         if (Math.abs(yInitial) <= functionValueAccuracy) {
  113.             return initial;
  114.         }

  116.         // Return the first endpoint if it is good enough.
  117.         final double yMin = func.applyAsDouble(min);
  118.         if (Math.abs(yMin) <= functionValueAccuracy) {
  119.             return min;
  120.         }

  122.         // Reduce interval if min and initial bracket the root.
  123.         if (Double.compare(yInitial * yMin, 0.0) < 0) {
  124.             return brent(func, min, initial, yMin, yInitial);
  125.         }

  127.         // Return the second endpoint if it is good enough.
  128.         final double yMax = func.applyAsDouble(max);
  129.         if (Math.abs(yMax) <= functionValueAccuracy) {
  130.             return max;
  131.         }

  133.         // Reduce interval if initial and max bracket the root.
  134.         if (Double.compare(yInitial * yMax, 0.0) < 0) {
  135.             return brent(func, initial, max, yInitial, yMax);
  136.         }

  138.         throw new SolverException(SolverException.BRACKETING, min, yMin, max, yMax);
  139.     }

  141.     /**
  142.      * Search for a zero inside the provided interval.
  143.      * This implementation is based on the algorithm described at page 58 of
  144.      * the book
  145.      * <blockquote>
  146.      *  <b>Algorithms for Minimization Without Derivatives</b>,
  147.      *  <i>Richard P. Brent</i>,
  148.      *  Dover 0-486-41998-3
  149.      * </blockquote>
  150.      *
  151.      * @param func Function to solve.
  152.      * @param lo Lower bound of the search interval.
  153.      * @param hi Higher bound of the search interval.
  154.      * @param fLo Function value at the lower bound of the search interval.
  155.      * @param fHi Function value at the higher bound of the search interval.
  156.      * @return the value where the function is zero.
  157.      */
  158.     private double brent(DoubleUnaryOperator func,
  159.                          double lo, double hi,
  160.                          double fLo, double fHi) {
  161.         double a = lo;
  162.         double fa = fLo;
  163.         double b = hi;
  164.         double fb = fHi;
  165.         double c = a;
  166.         double fc = fa;
  167.         double d = b - a;
  168.         double e = d;

  170.         final double t = absoluteAccuracy;
  171.         final double eps = relativeAccuracy;

  173.         while (true) {
  174.             if (Math.abs(fc) < Math.abs(fb)) {
  175.                 a = b;
  176.                 b = c;
  177.                 c = a;
  178.                 fa = fb;
  179.                 fb = fc;
  180.                 fc = fa;
  181.             }

  183.             final double tol = 2 * eps * Math.abs(b) + t;
  184.             final double m = 0.5 * (c - b);

  186.             if (Math.abs(m) <= tol ||
  187.                 equalsZero(fb))  {
  188.                 return b;
  189.             }
  190.             if (Math.abs(e) < tol ||
  191.                 Math.abs(fa) <= Math.abs(fb)) {
  192.                 // Force bisection.
  193.                 d = m;
  194.                 e = d;
  195.             } else {
  196.                 final double s = fb / fa;
  197.                 double p;
  198.                 double q;
  199.                 // The equality test (a == c) is intentional,
  200.                 // it is part of the original Brent's method and
  201.                 // it should NOT be replaced by proximity test.
  202.                 if (a == c) {
  203.                     // Linear interpolation.
  204.                     p = 2 * m * s;
  205.                     q = 1 - s;
  206.                 } else {
  207.                     // Inverse quadratic interpolation.
  208.                     q = fa / fc;
  209.                     final double r = fb / fc;
  210.                     p = s * (2 * m * q * (q - r) - (b - a) * (r - 1));
  211.                     q = (q - 1) * (r - 1) * (s - 1);
  212.                 }
  213.                 if (p > 0) {
  214.                     q = -q;
  215.                 } else {
  216.                     p = -p;
  217.                 }
  218.                 if (p >= 1.5 * m * q - Math.abs(tol * q) ||
  219.                     p >= Math.abs(0.5 * e * q)) {
  220.                     // Inverse quadratic interpolation gives a value
  221.                     // in the wrong direction, or progress is slow.
  222.                     // Fall back to bisection.
  223.                     d = m;
  224.                     e = d;
  225.                 } else {
  226.                     e = d;
  227.                     d = p / q;
  228.                 }
  229.             }
  230.             a = b;
  231.             fa = fb;

  233.             if (Math.abs(d) > tol) {
  234.                 b += d;
  235.             } else if (m > 0) {
  236.                 b += tol;
  237.             } else {
  238.                 b -= tol;
  239.             }
  240.             fb = func.applyAsDouble(b);
  241.             if ((fb > 0 && fc > 0) ||
  242.                 (fb <= 0 && fc <= 0)) {
  243.                 c = a;
  244.                 fc = fa;
  245.                 d = b - a;
  246.                 e = d;
  247.             }
  248.         }
  249.     }

  251.     /**
  252.      * Return true if the value is within 1 ULP of zero.
  253.      *
  254.      * @param value Value
  255.      * @return true if zero within a 1 ULP tolerance
  256.      */
  257.     private static boolean equalsZero(double value) {
  258.         return Math.abs(value) <= Double.MIN_VALUE;
  259.     }
  260. }
Created with JaCoCo 0.8.12.202403310830Code Coverage Report for Apache Commons Numbers 1.2