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.  *      https://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.  */

  18. /*
  19.  * This is not the original file distributed by the Apache Software Foundation
  20.  * It has been modified by the Hipparchus project
  21.  */
  22. package org.hipparchus.analysis.solvers;


  25. import org.hipparchus.exception.LocalizedCoreFormats;
  26. import org.hipparchus.exception.MathIllegalArgumentException;
  27. import org.hipparchus.exception.MathIllegalStateException;
  28. import org.hipparchus.util.FastMath;
  29. import org.hipparchus.util.Precision;

  31. /**
  32.  * This class implements the <a href="http://mathworld.wolfram.com/BrentsMethod.html">
  33.  * Brent algorithm</a> for finding zeros of real univariate functions.
  34.  * The function should be continuous but not necessarily smooth.
  35.  * The {@code solve} method returns a zero {@code x} of the function {@code f}
  36.  * in the given interval {@code [a, b]} to within a tolerance
  37.  * {@code 2 eps abs(x) + t} where {@code eps} is the relative accuracy and
  38.  * {@code t} is the absolute accuracy.
  39.  * <p>The given interval must bracket the root.</p>
  40.  * <p>
  41.  *  The reference implementation is given in chapter 4 of
  42.  *  <blockquote>
  43.  *   <b>Algorithms for Minimization Without Derivatives</b>,
  44.  *   <em>Richard P. Brent</em>,
  45.  *   Dover, 2002
  46.  *  </blockquote>
  47.  *
  48.  * @see BaseAbstractUnivariateSolver
  49.  */
  50. public class BrentSolver extends AbstractUnivariateSolver {

  52.     /** Default absolute accuracy. */
  53.     private static final double DEFAULT_ABSOLUTE_ACCURACY = 1e-6;

  55.     /**
  56.      * Construct a solver with default absolute accuracy (1e-6).
  57.      */
  58.     public BrentSolver() {
  59.         this(DEFAULT_ABSOLUTE_ACCURACY);
  60.     }
  61.     /**
  62.      * Construct a solver.
  63.      *
  64.      * @param absoluteAccuracy Absolute accuracy.
  65.      */
  66.     public BrentSolver(double absoluteAccuracy) {
  67.         super(absoluteAccuracy);
  68.     }
  69.     /**
  70.      * Construct a solver.
  71.      *
  72.      * @param relativeAccuracy Relative accuracy.
  73.      * @param absoluteAccuracy Absolute accuracy.
  74.      */
  75.     public BrentSolver(double relativeAccuracy,
  76.                        double absoluteAccuracy) {
  77.         super(relativeAccuracy, absoluteAccuracy);
  78.     }
  79.     /**
  80.      * Construct a solver.
  81.      *
  82.      * @param relativeAccuracy Relative accuracy.
  83.      * @param absoluteAccuracy Absolute accuracy.
  84.      * @param functionValueAccuracy Function value accuracy.
  85.      *
  86.      * @see BaseAbstractUnivariateSolver#BaseAbstractUnivariateSolver(double,double,double)
  87.      */
  88.     public BrentSolver(double relativeAccuracy,
  89.                        double absoluteAccuracy,
  90.                        double functionValueAccuracy) {
  91.         super(relativeAccuracy, absoluteAccuracy, functionValueAccuracy);
  92.     }

  94.     /**
  95.      * {@inheritDoc}
  96.      */
  97.     @Override
  98.     protected double doSolve()
  99.         throws MathIllegalArgumentException, MathIllegalStateException {
  100.         double min = getMin();
  101.         double max = getMax();
  102.         final double initial = getStartValue();
  103.         final double functionValueAccuracy = getFunctionValueAccuracy();

  105.         verifySequence(min, initial, max);

  107.         // Return the initial guess if it is good enough.
  108.         double yInitial = computeObjectiveValue(initial);
  109.         if (FastMath.abs(yInitial) <= functionValueAccuracy) {
  110.             return initial;
  111.         }

  113.         // Return the first endpoint if it is good enough.
  114.         double yMin = computeObjectiveValue(min);
  115.         if (FastMath.abs(yMin) <= functionValueAccuracy) {
  116.             return min;
  117.         }

  119.         // Reduce interval if min and initial bracket the root.
  120.         if (yInitial * yMin < 0) {
  121.             return brent(min, initial, yMin, yInitial);
  122.         }

  124.         // Return the second endpoint if it is good enough.
  125.         double yMax = computeObjectiveValue(max);
  126.         if (FastMath.abs(yMax) <= functionValueAccuracy) {
  127.             return max;
  128.         }

  130.         // Reduce interval if initial and max bracket the root.
  131.         if (yInitial * yMax < 0) {
  132.             return brent(initial, max, yInitial, yMax);
  133.         }

  135.         throw new MathIllegalArgumentException(LocalizedCoreFormats.NOT_BRACKETING_INTERVAL,
  136.                                                min, max, yMin, yMax);
  137.     }

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

  166.         final double t = getAbsoluteAccuracy();
  167.         final double eps = getRelativeAccuracy();

  169.         while (true) {
  170.             if (FastMath.abs(fc) < FastMath.abs(fb)) {
  171.                 a = b;
  172.                 b = c;
  173.                 c = a;
  174.                 fa = fb;
  175.                 fb = fc;
  176.                 fc = fa;
  177.             }

  179.             final double tol = 2 * eps * FastMath.abs(b) + t;
  180.             final double m = 0.5 * (c - b);

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

  230.             if (FastMath.abs(d) > tol) {
  231.                 b += d;
  232.             } else if (m > 0) {
  233.                 b += tol;
  234.             } else {
  235.                 b -= tol;
  236.             }
  237.             fb = computeObjectiveValue(b);
  238.             if ((fb > 0 && fc > 0) ||
  239.                 (fb <= 0 && fc <= 0)) {
  240.                 c = a;
  241.                 fc = fa;
  242.                 d = b - a;
  243.                 e = d;
  244.             }
  245.         }
  246.     }
  247. }
Created with JaCoCo 0.8.12.202403310830