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