001    /*
002     * Licensed to the Apache Software Foundation (ASF) under one or more
003     * contributor license agreements.  See the NOTICE file distributed with
004     * this work for additional information regarding copyright ownership.
005     * The ASF licenses this file to You under the Apache License, Version 2.0
006     * (the "License"); you may not use this file except in compliance with
007     * the License.  You may obtain a copy of the License at
008     *
009     *      http://www.apache.org/licenses/LICENSE-2.0
010     *
011     * Unless required by applicable law or agreed to in writing, software
012     * distributed under the License is distributed on an "AS IS" BASIS,
013     * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
014     * See the License for the specific language governing permissions and
015     * limitations under the License.
016     */
017    package org.apache.commons.math.analysis.solvers;
018    
019    
020    import org.apache.commons.math.FunctionEvaluationException;
021    import org.apache.commons.math.MathRuntimeException;
022    import org.apache.commons.math.MaxIterationsExceededException;
023    import org.apache.commons.math.analysis.UnivariateRealFunction;
024    import org.apache.commons.math.exception.util.LocalizedFormats;
025    import org.apache.commons.math.util.FastMath;
026    
027    /**
028     * Implements the <a href="http://mathworld.wolfram.com/BrentsMethod.html">
029     * Brent algorithm</a> for  finding zeros of real univariate functions.
030     * <p>
031     * The function should be continuous but not necessarily smooth.</p>
032     *
033     * @version $Revision:670469 $ $Date:2008-06-23 10:01:38 +0200 (lun., 23 juin 2008) $
034     */
035    public class BrentSolver extends UnivariateRealSolverImpl {
036    
037        /**
038         * Default absolute accuracy
039         * @since 2.1
040         */
041        public static final double DEFAULT_ABSOLUTE_ACCURACY = 1E-6;
042    
043        /** Default maximum number of iterations
044         * @since 2.1
045         */
046        public static final int DEFAULT_MAXIMUM_ITERATIONS = 100;
047    
048        /** Serializable version identifier */
049        private static final long serialVersionUID = 7694577816772532779L;
050    
051        /**
052         * Construct a solver for the given function.
053         *
054         * @param f function to solve.
055         * @deprecated as of 2.0 the function to solve is passed as an argument
056         * to the {@link #solve(UnivariateRealFunction, double, double)} or
057         * {@link UnivariateRealSolverImpl#solve(UnivariateRealFunction, double, double, double)}
058         * method.
059         */
060        @Deprecated
061        public BrentSolver(UnivariateRealFunction f) {
062            super(f, DEFAULT_MAXIMUM_ITERATIONS, DEFAULT_ABSOLUTE_ACCURACY);
063        }
064    
065        /**
066         * Construct a solver with default properties.
067         * @deprecated in 2.2 (to be removed in 3.0).
068         */
069        @Deprecated
070        public BrentSolver() {
071            super(DEFAULT_MAXIMUM_ITERATIONS, DEFAULT_ABSOLUTE_ACCURACY);
072        }
073    
074        /**
075         * Construct a solver with the given absolute accuracy.
076         *
077         * @param absoluteAccuracy lower bound for absolute accuracy of solutions returned by the solver
078         * @since 2.1
079         */
080        public BrentSolver(double absoluteAccuracy) {
081            super(DEFAULT_MAXIMUM_ITERATIONS, absoluteAccuracy);
082        }
083    
084        /**
085         * Contstruct a solver with the given maximum iterations and absolute accuracy.
086         *
087         * @param maximumIterations maximum number of iterations
088         * @param absoluteAccuracy lower bound for absolute accuracy of solutions returned by the solver
089         * @since 2.1
090         */
091        public BrentSolver(int maximumIterations, double absoluteAccuracy) {
092            super(maximumIterations, absoluteAccuracy);
093        }
094    
095        /** {@inheritDoc} */
096        @Deprecated
097        public double solve(double min, double max)
098            throws MaxIterationsExceededException, FunctionEvaluationException {
099            return solve(f, min, max);
100        }
101    
102        /** {@inheritDoc} */
103        @Deprecated
104        public double solve(double min, double max, double initial)
105            throws MaxIterationsExceededException, FunctionEvaluationException {
106            return solve(f, min, max, initial);
107        }
108    
109        /**
110         * Find a zero in the given interval with an initial guess.
111         * <p>Throws <code>IllegalArgumentException</code> if the values of the
112         * function at the three points have the same sign (note that it is
113         * allowed to have endpoints with the same sign if the initial point has
114         * opposite sign function-wise).</p>
115         *
116         * @param f function to solve.
117         * @param min the lower bound for the interval.
118         * @param max the upper bound for the interval.
119         * @param initial the start value to use (must be set to min if no
120         * initial point is known).
121         * @return the value where the function is zero
122         * @throws MaxIterationsExceededException the maximum iteration count is exceeded
123         * @throws FunctionEvaluationException if an error occurs evaluating  the function
124         * @throws IllegalArgumentException if initial is not between min and max
125         * (even if it <em>is</em> a root)
126         * @deprecated in 2.2 (to be removed in 3.0).
127         */
128        @Deprecated
129        public double solve(final UnivariateRealFunction f,
130                            final double min, final double max, final double initial)
131            throws MaxIterationsExceededException, FunctionEvaluationException {
132    
133            clearResult();
134            if ((initial < min) || (initial > max)) {
135                throw MathRuntimeException.createIllegalArgumentException(
136                      LocalizedFormats.INVALID_INTERVAL_INITIAL_VALUE_PARAMETERS,
137                      min, initial, max);
138            }
139    
140            // return the initial guess if it is good enough
141            double yInitial = f.value(initial);
142            if (FastMath.abs(yInitial) <= functionValueAccuracy) {
143                setResult(initial, 0);
144                return result;
145            }
146    
147            // return the first endpoint if it is good enough
148            double yMin = f.value(min);
149            if (FastMath.abs(yMin) <= functionValueAccuracy) {
150                setResult(min, 0);
151                return result;
152            }
153    
154            // reduce interval if min and initial bracket the root
155            if (yInitial * yMin < 0) {
156                return solve(f, min, yMin, initial, yInitial, min, yMin);
157            }
158    
159            // return the second endpoint if it is good enough
160            double yMax = f.value(max);
161            if (FastMath.abs(yMax) <= functionValueAccuracy) {
162                setResult(max, 0);
163                return result;
164            }
165    
166            // reduce interval if initial and max bracket the root
167            if (yInitial * yMax < 0) {
168                return solve(f, initial, yInitial, max, yMax, initial, yInitial);
169            }
170    
171            throw MathRuntimeException.createIllegalArgumentException(
172                  LocalizedFormats.SAME_SIGN_AT_ENDPOINTS, min, max, yMin, yMax);
173    
174        }
175    
176        /**
177         * Find a zero in the given interval with an initial guess.
178         * <p>Throws <code>IllegalArgumentException</code> if the values of the
179         * function at the three points have the same sign (note that it is
180         * allowed to have endpoints with the same sign if the initial point has
181         * opposite sign function-wise).</p>
182         *
183         * @param f function to solve.
184         * @param min the lower bound for the interval.
185         * @param max the upper bound for the interval.
186         * @param initial the start value to use (must be set to min if no
187         * initial point is known).
188         * @param maxEval Maximum number of evaluations.
189         * @return the value where the function is zero
190         * @throws MaxIterationsExceededException the maximum iteration count is exceeded
191         * @throws FunctionEvaluationException if an error occurs evaluating  the function
192         * @throws IllegalArgumentException if initial is not between min and max
193         * (even if it <em>is</em> a root)
194         */
195        @Override
196        public double solve(int maxEval, final UnivariateRealFunction f,
197                            final double min, final double max, final double initial)
198            throws MaxIterationsExceededException, FunctionEvaluationException {
199            setMaximalIterationCount(maxEval);
200            return solve(f, min, max, initial);
201        }
202    
203        /**
204         * Find a zero in the given interval.
205         * <p>
206         * Requires that the values of the function at the endpoints have opposite
207         * signs. An <code>IllegalArgumentException</code> is thrown if this is not
208         * the case.</p>
209         *
210         * @param f the function to solve
211         * @param min the lower bound for the interval.
212         * @param max the upper bound for the interval.
213         * @return the value where the function is zero
214         * @throws MaxIterationsExceededException if the maximum iteration count is exceeded
215         * @throws FunctionEvaluationException if an error occurs evaluating the function
216         * @throws IllegalArgumentException if min is not less than max or the
217         * signs of the values of the function at the endpoints are not opposites
218         * @deprecated in 2.2 (to be removed in 3.0).
219         */
220        @Deprecated
221        public double solve(final UnivariateRealFunction f,
222                            final double min, final double max)
223            throws MaxIterationsExceededException, FunctionEvaluationException {
224    
225            clearResult();
226            verifyInterval(min, max);
227    
228            double ret = Double.NaN;
229    
230            double yMin = f.value(min);
231            double yMax = f.value(max);
232    
233            // Verify bracketing
234            double sign = yMin * yMax;
235            if (sign > 0) {
236                // check if either value is close to a zero
237                if (FastMath.abs(yMin) <= functionValueAccuracy) {
238                    setResult(min, 0);
239                    ret = min;
240                } else if (FastMath.abs(yMax) <= functionValueAccuracy) {
241                    setResult(max, 0);
242                    ret = max;
243                } else {
244                    // neither value is close to zero and min and max do not bracket root.
245                    throw MathRuntimeException.createIllegalArgumentException(
246                            LocalizedFormats.SAME_SIGN_AT_ENDPOINTS, min, max, yMin, yMax);
247                }
248            } else if (sign < 0){
249                // solve using only the first endpoint as initial guess
250                ret = solve(f, min, yMin, max, yMax, min, yMin);
251            } else {
252                // either min or max is a root
253                if (yMin == 0.0) {
254                    ret = min;
255                } else {
256                    ret = max;
257                }
258            }
259    
260            return ret;
261        }
262    
263        /**
264         * Find a zero in the given interval.
265         * <p>
266         * Requires that the values of the function at the endpoints have opposite
267         * signs. An <code>IllegalArgumentException</code> is thrown if this is not
268         * the case.</p>
269         *
270         * @param f the function to solve
271         * @param min the lower bound for the interval.
272         * @param max the upper bound for the interval.
273         * @param maxEval Maximum number of evaluations.
274         * @return the value where the function is zero
275         * @throws MaxIterationsExceededException if the maximum iteration count is exceeded
276         * @throws FunctionEvaluationException if an error occurs evaluating the function
277         * @throws IllegalArgumentException if min is not less than max or the
278         * signs of the values of the function at the endpoints are not opposites
279         */
280        @Override
281        public double solve(int maxEval, final UnivariateRealFunction f,
282                            final double min, final double max)
283            throws MaxIterationsExceededException, FunctionEvaluationException {
284            setMaximalIterationCount(maxEval);
285            return solve(f, min, max);
286        }
287    
288        /**
289         * Find a zero starting search according to the three provided points.
290         * @param f the function to solve
291         * @param x0 old approximation for the root
292         * @param y0 function value at the approximation for the root
293         * @param x1 last calculated approximation for the root
294         * @param y1 function value at the last calculated approximation
295         * for the root
296         * @param x2 bracket point (must be set to x0 if no bracket point is
297         * known, this will force starting with linear interpolation)
298         * @param y2 function value at the bracket point.
299         * @return the value where the function is zero
300         * @throws MaxIterationsExceededException if the maximum iteration count is exceeded
301         * @throws FunctionEvaluationException if an error occurs evaluating the function
302         */
303        private double solve(final UnivariateRealFunction f,
304                             double x0, double y0,
305                             double x1, double y1,
306                             double x2, double y2)
307        throws MaxIterationsExceededException, FunctionEvaluationException {
308    
309            double delta = x1 - x0;
310            double oldDelta = delta;
311    
312            int i = 0;
313            while (i < maximalIterationCount) {
314                if (FastMath.abs(y2) < FastMath.abs(y1)) {
315                    // use the bracket point if is better than last approximation
316                    x0 = x1;
317                    x1 = x2;
318                    x2 = x0;
319                    y0 = y1;
320                    y1 = y2;
321                    y2 = y0;
322                }
323                if (FastMath.abs(y1) <= functionValueAccuracy) {
324                    // Avoid division by very small values. Assume
325                    // the iteration has converged (the problem may
326                    // still be ill conditioned)
327                    setResult(x1, i);
328                    return result;
329                }
330                double dx = x2 - x1;
331                double tolerance =
332                    FastMath.max(relativeAccuracy * FastMath.abs(x1), absoluteAccuracy);
333                if (FastMath.abs(dx) <= tolerance) {
334                    setResult(x1, i);
335                    return result;
336                }
337                if ((FastMath.abs(oldDelta) < tolerance) ||
338                        (FastMath.abs(y0) <= FastMath.abs(y1))) {
339                    // Force bisection.
340                    delta = 0.5 * dx;
341                    oldDelta = delta;
342                } else {
343                    double r3 = y1 / y0;
344                    double p;
345                    double p1;
346                    // the equality test (x0 == x2) is intentional,
347                    // it is part of the original Brent's method,
348                    // it should NOT be replaced by proximity test
349                    if (x0 == x2) {
350                        // Linear interpolation.
351                        p = dx * r3;
352                        p1 = 1.0 - r3;
353                    } else {
354                        // Inverse quadratic interpolation.
355                        double r1 = y0 / y2;
356                        double r2 = y1 / y2;
357                        p = r3 * (dx * r1 * (r1 - r2) - (x1 - x0) * (r2 - 1.0));
358                        p1 = (r1 - 1.0) * (r2 - 1.0) * (r3 - 1.0);
359                    }
360                    if (p > 0.0) {
361                        p1 = -p1;
362                    } else {
363                        p = -p;
364                    }
365                    if (2.0 * p >= 1.5 * dx * p1 - FastMath.abs(tolerance * p1) ||
366                            p >= FastMath.abs(0.5 * oldDelta * p1)) {
367                        // Inverse quadratic interpolation gives a value
368                        // in the wrong direction, or progress is slow.
369                        // Fall back to bisection.
370                        delta = 0.5 * dx;
371                        oldDelta = delta;
372                    } else {
373                        oldDelta = delta;
374                        delta = p / p1;
375                    }
376                }
377                // Save old X1, Y1
378                x0 = x1;
379                y0 = y1;
380                // Compute new X1, Y1
381                if (FastMath.abs(delta) > tolerance) {
382                    x1 = x1 + delta;
383                } else if (dx > 0.0) {
384                    x1 = x1 + 0.5 * tolerance;
385                } else if (dx <= 0.0) {
386                    x1 = x1 - 0.5 * tolerance;
387                }
388                y1 = f.value(x1);
389                if ((y1 > 0) == (y2 > 0)) {
390                    x2 = x0;
391                    y2 = y0;
392                    delta = x1 - x0;
393                    oldDelta = delta;
394                }
395                i++;
396            }
397            throw new MaxIterationsExceededException(maximalIterationCount);
398        }
399    }