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