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