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    *      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.math3.analysis.solvers;
18  
19  
20  import org.apache.commons.math3.exception.NoBracketingException;
21  import org.apache.commons.math3.exception.TooManyEvaluationsException;
22  import org.apache.commons.math3.exception.NumberIsTooLargeException;
23  import org.apache.commons.math3.util.FastMath;
24  import org.apache.commons.math3.util.Precision;
25  
26  /**
27   * This class implements the <a href="http://mathworld.wolfram.com/BrentsMethod.html">
28   * Brent algorithm</a> for finding zeros of real univariate functions.
29   * The function should be continuous but not necessarily smooth.
30   * The {@code solve} method returns a zero {@code x} of the function {@code f}
31   * in the given interval {@code [a, b]} to within a tolerance
32   * {@code 6 eps abs(x) + t} where {@code eps} is the relative accuracy and
33   * {@code t} is the absolute accuracy.
34   * The given interval must bracket the root.
35   *
36   * @version $Id: BrentSolver.java 1379560 2012-08-31 19:40:30Z erans $
37   */
38  public class BrentSolver extends AbstractUnivariateSolver {
39  
40      /** Default absolute accuracy. */
41      private static final double DEFAULT_ABSOLUTE_ACCURACY = 1e-6;
42  
43      /**
44       * Construct a solver with default accuracy (1e-6).
45       */
46      public BrentSolver() {
47          this(DEFAULT_ABSOLUTE_ACCURACY);
48      }
49      /**
50       * Construct a solver.
51       *
52       * @param absoluteAccuracy Absolute accuracy.
53       */
54      public BrentSolver(double absoluteAccuracy) {
55          super(absoluteAccuracy);
56      }
57      /**
58       * Construct a solver.
59       *
60       * @param relativeAccuracy Relative accuracy.
61       * @param absoluteAccuracy Absolute accuracy.
62       */
63      public BrentSolver(double relativeAccuracy,
64                         double absoluteAccuracy) {
65          super(relativeAccuracy, absoluteAccuracy);
66      }
67      /**
68       * Construct a solver.
69       *
70       * @param relativeAccuracy Relative accuracy.
71       * @param absoluteAccuracy Absolute accuracy.
72       * @param functionValueAccuracy Function value accuracy.
73       */
74      public BrentSolver(double relativeAccuracy,
75                         double absoluteAccuracy,
76                         double functionValueAccuracy) {
77          super(relativeAccuracy, absoluteAccuracy, functionValueAccuracy);
78      }
79  
80      /**
81       * {@inheritDoc}
82       */
83      @Override
84      protected double doSolve()
85          throws NoBracketingException,
86                 TooManyEvaluationsException,
87                 NumberIsTooLargeException {
88          double min = getMin();
89          double max = getMax();
90          final double initial = getStartValue();
91          final double functionValueAccuracy = getFunctionValueAccuracy();
92  
93          verifySequence(min, initial, max);
94  
95          // Return the initial guess if it is good enough.
96          double yInitial = computeObjectiveValue(initial);
97          if (FastMath.abs(yInitial) <= functionValueAccuracy) {
98              return initial;
99          }
100 
101         // Return the first endpoint if it is good enough.
102         double yMin = computeObjectiveValue(min);
103         if (FastMath.abs(yMin) <= functionValueAccuracy) {
104             return min;
105         }
106 
107         // Reduce interval if min and initial bracket the root.
108         if (yInitial * yMin < 0) {
109             return brent(min, initial, yMin, yInitial);
110         }
111 
112         // Return the second endpoint if it is good enough.
113         double yMax = computeObjectiveValue(max);
114         if (FastMath.abs(yMax) <= functionValueAccuracy) {
115             return max;
116         }
117 
118         // Reduce interval if initial and max bracket the root.
119         if (yInitial * yMax < 0) {
120             return brent(initial, max, yInitial, yMax);
121         }
122 
123         throw new NoBracketingException(min, max, yMin, yMax);
124     }
125 
126     /**
127      * Search for a zero inside the provided interval.
128      * This implementation is based on the algorithm described at page 58 of
129      * the book
130      * <quote>
131      *  <b>Algorithms for Minimization Without Derivatives</b>
132      *  <it>Richard P. Brent</it>
133      *  Dover 0-486-41998-3
134      * </quote>
135      *
136      * @param lo Lower bound of the search interval.
137      * @param hi Higher bound of the search interval.
138      * @param fLo Function value at the lower bound of the search interval.
139      * @param fHi Function value at the higher bound of the search interval.
140      * @return the value where the function is zero.
141      */
142     private double brent(double lo, double hi,
143                          double fLo, double fHi) {
144         double a = lo;
145         double fa = fLo;
146         double b = hi;
147         double fb = fHi;
148         double c = a;
149         double fc = fa;
150         double d = b - a;
151         double e = d;
152 
153         final double t = getAbsoluteAccuracy();
154         final double eps = getRelativeAccuracy();
155 
156         while (true) {
157             if (FastMath.abs(fc) < FastMath.abs(fb)) {
158                 a = b;
159                 b = c;
160                 c = a;
161                 fa = fb;
162                 fb = fc;
163                 fc = fa;
164             }
165 
166             final double tol = 2 * eps * FastMath.abs(b) + t;
167             final double m = 0.5 * (c - b);
168 
169             if (FastMath.abs(m) <= tol ||
170                 Precision.equals(fb, 0))  {
171                 return b;
172             }
173             if (FastMath.abs(e) < tol ||
174                 FastMath.abs(fa) <= FastMath.abs(fb)) {
175                 // Force bisection.
176                 d = m;
177                 e = d;
178             } else {
179                 double s = fb / fa;
180                 double p;
181                 double q;
182                 // The equality test (a == c) is intentional,
183                 // it is part of the original Brent's method and
184                 // it should NOT be replaced by proximity test.
185                 if (a == c) {
186                     // Linear interpolation.
187                     p = 2 * m * s;
188                     q = 1 - s;
189                 } else {
190                     // Inverse quadratic interpolation.
191                     q = fa / fc;
192                     final double r = fb / fc;
193                     p = s * (2 * m * q * (q - r) - (b - a) * (r - 1));
194                     q = (q - 1) * (r - 1) * (s - 1);
195                 }
196                 if (p > 0) {
197                     q = -q;
198                 } else {
199                     p = -p;
200                 }
201                 s = e;
202                 e = d;
203                 if (p >= 1.5 * m * q - FastMath.abs(tol * q) ||
204                     p >= FastMath.abs(0.5 * s * q)) {
205                     // Inverse quadratic interpolation gives a value
206                     // in the wrong direction, or progress is slow.
207                     // Fall back to bisection.
208                     d = m;
209                     e = d;
210                 } else {
211                     d = p / q;
212                 }
213             }
214             a = b;
215             fa = fb;
216 
217             if (FastMath.abs(d) > tol) {
218                 b += d;
219             } else if (m > 0) {
220                 b += tol;
221             } else {
222                 b -= tol;
223             }
224             fb = computeObjectiveValue(b);
225             if ((fb > 0 && fc > 0) ||
226                 (fb <= 0 && fc <= 0)) {
227                 c = a;
228                 fc = fa;
229                 d = b - a;
230                 e = d;
231             }
232         }
233     }
234 }