View Javadoc
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   */
17  
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;
23  
24  
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;
30  
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 {
51  
52      /** Default absolute accuracy. */
53      private static final double DEFAULT_ABSOLUTE_ACCURACY = 1e-6;
54  
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      }
93  
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();
104 
105         verifySequence(min, initial, max);
106 
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         }
112 
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         }
118 
119         // Reduce interval if min and initial bracket the root.
120         if (yInitial * yMin < 0) {
121             return brent(min, initial, yMin, yInitial);
122         }
123 
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         }
129 
130         // Reduce interval if initial and max bracket the root.
131         if (yInitial * yMax < 0) {
132             return brent(initial, max, yInitial, yMax);
133         }
134 
135         throw new MathIllegalArgumentException(LocalizedCoreFormats.NOT_BRACKETING_INTERVAL,
136                                                min, max, yMin, yMax);
137     }
138 
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;
165 
166         final double t = getAbsoluteAccuracy();
167         final double eps = getRelativeAccuracy();
168 
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             }
178 
179             final double tol = 2 * eps * FastMath.abs(b) + t;
180             final double m = 0.5 * (c - b);
181 
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;
229 
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 }