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.math.analysis;
18  
19  
20  import org.apache.commons.math.FunctionEvaluationException;
21  import org.apache.commons.math.MaxIterationsExceededException;
22  
23  /**
24   * Implements the <a href="http://mathworld.wolfram.com/BrentsMethod.html">
25   * Brent algorithm</a> for  finding zeros of real univariate functions.
26   * <p>
27   * The function should be continuous but not necessarily smooth.</p>
28   *  
29   * @version $Revision: 617826 $ $Date: 2008-02-02 09:39:18 -0700 (Sat, 02 Feb 2008) $
30   */
31  public class BrentSolver extends UnivariateRealSolverImpl {
32      
33      /** Serializable version identifier */
34      private static final long serialVersionUID = -2136672307739067002L;
35  
36      /**
37       * Construct a solver for the given function.
38       * 
39       * @param f function to solve.
40       */
41      public BrentSolver(UnivariateRealFunction f) {
42          super(f, 100, 1E-6);
43      }
44  
45      /**
46       * Find a zero in the given interval with an initial guess.
47       * <p>Throws <code>IllegalArgumentException</code> if the values of the
48       * function at the three points have the same sign (note that it is
49       * allowed to have endpoints with the same sign if the initial point has
50       * opposite sign function-wise).</p>
51       * 
52       * @param min the lower bound for the interval.
53       * @param max the upper bound for the interval.
54       * @param initial the start value to use (must be set to min if no
55       * initial point is known).
56       * @return the value where the function is zero
57       * @throws MaxIterationsExceededException the maximum iteration count
58       * is exceeded 
59       * @throws FunctionEvaluationException if an error occurs evaluating
60       *  the function
61       * @throws IllegalArgumentException if initial is not between min and max
62       * (even if it <em>is</em> a root)
63       */
64      public double solve(double min, double max, double initial)
65          throws MaxIterationsExceededException, FunctionEvaluationException {
66  
67          if (((initial - min) * (max -initial)) < 0) {
68              throw new IllegalArgumentException("Initial guess is not in search" +
69                        " interval." + "  Initial: " + initial +
70                        "  Endpoints: [" + min + "," + max + "]");
71          }
72  
73          // return the initial guess if it is good enough
74          double yInitial = f.value(initial);
75          if (Math.abs(yInitial) <= functionValueAccuracy) {
76              setResult(initial, 0);
77              return result;
78          }
79  
80          // return the first endpoint if it is good enough
81          double yMin = f.value(min);
82          if (Math.abs(yMin) <= functionValueAccuracy) {
83              setResult(yMin, 0);
84              return result;
85          }
86  
87          // reduce interval if min and initial bracket the root
88          if (yInitial * yMin < 0) {
89              return solve(min, yMin, initial, yInitial, min, yMin);
90          }
91  
92          // return the second endpoint if it is good enough
93          double yMax = f.value(max);
94          if (Math.abs(yMax) <= functionValueAccuracy) {
95              setResult(yMax, 0);
96              return result;
97          }
98  
99          // reduce interval if initial and max bracket the root
100         if (yInitial * yMax < 0) {
101             return solve(initial, yInitial, max, yMax, initial, yInitial);
102         }
103 
104         // full Brent algorithm starting with provided initial guess
105         return solve(min, yMin, max, yMax, initial, yInitial);
106 
107     }
108     
109     /**
110      * Find a zero in the given interval.
111      * <p>
112      * Requires that the values of the function at the endpoints have opposite
113      * signs. An <code>IllegalArgumentException</code> is thrown if this is not
114      * the case.</p>
115      * 
116      * @param min the lower bound for the interval.
117      * @param max the upper bound for the interval.
118      * @return the value where the function is zero
119      * @throws MaxIterationsExceededException if the maximum iteration count is exceeded
120      * @throws FunctionEvaluationException if an error occurs evaluating the
121      * function 
122      * @throws IllegalArgumentException if min is not less than max or the
123      * signs of the values of the function at the endpoints are not opposites
124      */
125     public double solve(double min, double max) throws MaxIterationsExceededException, 
126         FunctionEvaluationException {
127         
128         clearResult();
129         verifyInterval(min, max);
130         
131         double yMin = f.value(min);
132         double yMax = f.value(max);
133         
134         // Verify bracketing
135         if (yMin * yMax >= 0) {
136             throw new IllegalArgumentException
137             ("Function values at endpoints do not have different signs." +
138                     "  Endpoints: [" + min + "," + max + "]" + 
139                     "  Values: [" + yMin + "," + yMax + "]");       
140         }
141 
142         // solve using only the first endpoint as initial guess
143         return solve(min, yMin, max, yMax, min, yMin);
144 
145     }
146         
147     /**
148      * Find a zero starting search according to the three provided points.
149      * @param x0 old approximation for the root
150      * @param y0 function value at the approximation for the root
151      * @param x1 last calculated approximation for the root
152      * @param y1 function value at the last calculated approximation
153      * for the root
154      * @param x2 bracket point (must be set to x0 if no bracket point is
155      * known, this will force starting with linear interpolation)
156      * @param y2 function value at the bracket point.
157      * @return the value where the function is zero
158      * @throws MaxIterationsExceededException if the maximum iteration count
159      * is exceeded
160      * @throws FunctionEvaluationException if an error occurs evaluating
161      * the function 
162      */
163     private double solve(double x0, double y0,
164                          double x1, double y1,
165                          double x2, double y2)
166     throws MaxIterationsExceededException, FunctionEvaluationException {
167 
168         double delta = x1 - x0;
169         double oldDelta = delta;
170 
171         int i = 0;
172         while (i < maximalIterationCount) {
173             if (Math.abs(y2) < Math.abs(y1)) {
174                 // use the bracket point if is better than last approximation
175                 x0 = x1;
176                 x1 = x2;
177                 x2 = x0;
178                 y0 = y1;
179                 y1 = y2;
180                 y2 = y0;
181             }
182             if (Math.abs(y1) <= functionValueAccuracy) {
183                 // Avoid division by very small values. Assume
184                 // the iteration has converged (the problem may
185                 // still be ill conditioned)
186                 setResult(x1, i);
187                 return result;
188             }
189             double dx = (x2 - x1);
190             double tolerance =
191                 Math.max(relativeAccuracy * Math.abs(x1), absoluteAccuracy);
192             if (Math.abs(dx) <= tolerance) {
193                 setResult(x1, i);
194                 return result;
195             }
196             if ((Math.abs(oldDelta) < tolerance) ||
197                     (Math.abs(y0) <= Math.abs(y1))) {
198                 // Force bisection.
199                 delta = 0.5 * dx;
200                 oldDelta = delta;
201             } else {
202                 double r3 = y1 / y0;
203                 double p;
204                 double p1;
205                 // the equality test (x0 == x2) is intentional,
206                 // it is part of the original Brent's method,
207                 // it should NOT be replaced by proximity test
208                 if (x0 == x2) {
209                     // Linear interpolation.
210                     p = dx * r3;
211                     p1 = 1.0 - r3;
212                 } else {
213                     // Inverse quadratic interpolation.
214                     double r1 = y0 / y2;
215                     double r2 = y1 / y2;
216                     p = r3 * (dx * r1 * (r1 - r2) - (x1 - x0) * (r2 - 1.0));
217                     p1 = (r1 - 1.0) * (r2 - 1.0) * (r3 - 1.0);
218                 }
219                 if (p > 0.0) {
220                     p1 = -p1;
221                 } else {
222                     p = -p;
223                 }
224                 if (2.0 * p >= 1.5 * dx * p1 - Math.abs(tolerance * p1) ||
225                         p >= Math.abs(0.5 * oldDelta * p1)) {
226                     // Inverse quadratic interpolation gives a value
227                     // in the wrong direction, or progress is slow.
228                     // Fall back to bisection.
229                     delta = 0.5 * dx;
230                     oldDelta = delta;
231                 } else {
232                     oldDelta = delta;
233                     delta = p / p1;
234                 }
235             }
236             // Save old X1, Y1 
237             x0 = x1;
238             y0 = y1;
239             // Compute new X1, Y1
240             if (Math.abs(delta) > tolerance) {
241                 x1 = x1 + delta;
242             } else if (dx > 0.0) {
243                 x1 = x1 + 0.5 * tolerance;
244             } else if (dx <= 0.0) {
245                 x1 = x1 - 0.5 * tolerance;
246             }
247             y1 = f.value(x1);
248             if ((y1 > 0) == (y2 > 0)) {
249                 x2 = x0;
250                 y2 = y0;
251                 delta = x1 - x0;
252                 oldDelta = delta;
253             }
254             i++;
255         }
256         throw new MaxIterationsExceededException(maximalIterationCount);
257     }
258 }