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    
025    /**
026     * Implements the <a href="http://mathworld.wolfram.com/BrentsMethod.html">
027     * Brent algorithm</a> for  finding zeros of real univariate functions.
028     * <p>
029     * The function should be continuous but not necessarily smooth.</p>
030     *
031     * @version $Revision:670469 $ $Date:2008-06-23 10:01:38 +0200 (lun., 23 juin 2008) $
032     */
033    public class BrentSolver extends UnivariateRealSolverImpl {
034    
035        /**
036         * Default absolute accuracy
037         * @since 2.1
038         */
039        public static final double DEFAULT_ABSOLUTE_ACCURACY = 1E-6;
040    
041        /** Default maximum number of iterations
042         * @since 2.1
043         */
044        public static final int DEFAULT_MAXIMUM_ITERATIONS = 100;
045    
046        /** Error message for non-bracketing interval. */
047        private static final String NON_BRACKETING_MESSAGE =
048            "function values at endpoints do not have different signs.  " +
049            "Endpoints: [{0}, {1}], Values: [{2}, {3}]";
050    
051        /** Serializable version identifier */
052        private static final long serialVersionUID = 7694577816772532779L;
053    
054        /**
055         * Construct a solver for the given function.
056         *
057         * @param f function to solve.
058         * @deprecated as of 2.0 the function to solve is passed as an argument
059         * to the {@link #solve(UnivariateRealFunction, double, double)} or
060         * {@link UnivariateRealSolverImpl#solve(UnivariateRealFunction, double, double, double)}
061         * method.
062         */
063        @Deprecated
064        public BrentSolver(UnivariateRealFunction f) {
065            super(f, DEFAULT_MAXIMUM_ITERATIONS, DEFAULT_ABSOLUTE_ACCURACY);
066        }
067    
068        /**
069         * Construct a solver with default properties.
070         */
071        public BrentSolver() {
072            super(DEFAULT_MAXIMUM_ITERATIONS, DEFAULT_ABSOLUTE_ACCURACY);
073        }
074    
075        /**
076         * Construct a solver with the given absolute accuracy.
077         *
078         * @param absoluteAccuracy lower bound for absolute accuracy of solutions returned by the solver
079         * @since 2.1
080         */
081        public BrentSolver(double absoluteAccuracy) {
082            super(DEFAULT_MAXIMUM_ITERATIONS, absoluteAccuracy);
083        }
084    
085        /**
086         * Contstruct a solver with the given maximum iterations and absolute accuracy.
087         *
088         * @param maximumIterations maximum number of iterations
089         * @param absoluteAccuracy lower bound for absolute accuracy of solutions returned by the solver
090         * @since 2.1
091         */
092        public BrentSolver(int maximumIterations, double absoluteAccuracy) {
093            super(maximumIterations, absoluteAccuracy);
094        }
095    
096        /** {@inheritDoc} */
097        @Deprecated
098        public double solve(double min, double max)
099            throws MaxIterationsExceededException, FunctionEvaluationException {
100            return solve(f, min, max);
101        }
102    
103        /** {@inheritDoc} */
104        @Deprecated
105        public double solve(double min, double max, double initial)
106            throws MaxIterationsExceededException, FunctionEvaluationException {
107            return solve(f, min, max, initial);
108        }
109    
110        /**
111         * Find a zero in the given interval with an initial guess.
112         * <p>Throws <code>IllegalArgumentException</code> if the values of the
113         * function at the three points have the same sign (note that it is
114         * allowed to have endpoints with the same sign if the initial point has
115         * opposite sign function-wise).</p>
116         *
117         * @param f function to solve.
118         * @param min the lower bound for the interval.
119         * @param max the upper bound for the interval.
120         * @param initial the start value to use (must be set to min if no
121         * initial point is known).
122         * @return the value where the function is zero
123         * @throws MaxIterationsExceededException the maximum iteration count
124         * is exceeded
125         * @throws FunctionEvaluationException if an error occurs evaluating
126         *  the function
127         * @throws IllegalArgumentException if initial is not between min and max
128         * (even if it <em>is</em> a root)
129         */
130        public double solve(final UnivariateRealFunction f,
131                            final double min, final double max, final double initial)
132            throws MaxIterationsExceededException, FunctionEvaluationException {
133    
134            clearResult();
135            if ((initial < min) || (initial > max)) {
136                throw MathRuntimeException.createIllegalArgumentException(
137                      "invalid interval, initial value parameters:  lower={0}, initial={1}, upper={2}",
138                      min, initial, max);
139            }
140    
141            // return the initial guess if it is good enough
142            double yInitial = f.value(initial);
143            if (Math.abs(yInitial) <= functionValueAccuracy) {
144                setResult(initial, 0);
145                return result;
146            }
147    
148            // return the first endpoint if it is good enough
149            double yMin = f.value(min);
150            if (Math.abs(yMin) <= functionValueAccuracy) {
151                setResult(min, 0);
152                return result;
153            }
154    
155            // reduce interval if min and initial bracket the root
156            if (yInitial * yMin < 0) {
157                return solve(f, min, yMin, initial, yInitial, min, yMin);
158            }
159    
160            // return the second endpoint if it is good enough
161            double yMax = f.value(max);
162            if (Math.abs(yMax) <= functionValueAccuracy) {
163                setResult(max, 0);
164                return result;
165            }
166    
167            // reduce interval if initial and max bracket the root
168            if (yInitial * yMax < 0) {
169                return solve(f, initial, yInitial, max, yMax, initial, yInitial);
170            }
171    
172            throw MathRuntimeException.createIllegalArgumentException(
173                  NON_BRACKETING_MESSAGE, min, max, yMin, yMax);
174    
175        }
176    
177        /**
178         * Find a zero in the given interval.
179         * <p>
180         * Requires that the values of the function at the endpoints have opposite
181         * signs. An <code>IllegalArgumentException</code> is thrown if this is not
182         * the case.</p>
183         *
184         * @param f the function to solve
185         * @param min the lower bound for the interval.
186         * @param max the upper bound for the interval.
187         * @return the value where the function is zero
188         * @throws MaxIterationsExceededException if the maximum iteration count is exceeded
189         * @throws FunctionEvaluationException if an error occurs evaluating the
190         * function
191         * @throws IllegalArgumentException if min is not less than max or the
192         * signs of the values of the function at the endpoints are not opposites
193         */
194        public double solve(final UnivariateRealFunction f,
195                            final double min, final double max)
196            throws MaxIterationsExceededException,
197            FunctionEvaluationException {
198    
199            clearResult();
200            verifyInterval(min, max);
201    
202            double ret = Double.NaN;
203    
204            double yMin = f.value(min);
205            double yMax = f.value(max);
206    
207            // Verify bracketing
208            double sign = yMin * yMax;
209            if (sign > 0) {
210                // check if either value is close to a zero
211                if (Math.abs(yMin) <= functionValueAccuracy) {
212                    setResult(min, 0);
213                    ret = min;
214                } else if (Math.abs(yMax) <= functionValueAccuracy) {
215                    setResult(max, 0);
216                    ret = max;
217                } else {
218                    // neither value is close to zero and min and max do not bracket root.
219                    throw MathRuntimeException.createIllegalArgumentException(
220                            NON_BRACKETING_MESSAGE, min, max, yMin, yMax);
221                }
222            } else if (sign < 0){
223                // solve using only the first endpoint as initial guess
224                ret = solve(f, min, yMin, max, yMax, min, yMin);
225            } else {
226                // either min or max is a root
227                if (yMin == 0.0) {
228                    ret = min;
229                } else {
230                    ret = max;
231                }
232            }
233    
234            return ret;
235        }
236    
237        /**
238         * Find a zero starting search according to the three provided points.
239         * @param f the function to solve
240         * @param x0 old approximation for the root
241         * @param y0 function value at the approximation for the root
242         * @param x1 last calculated approximation for the root
243         * @param y1 function value at the last calculated approximation
244         * for the root
245         * @param x2 bracket point (must be set to x0 if no bracket point is
246         * known, this will force starting with linear interpolation)
247         * @param y2 function value at the bracket point.
248         * @return the value where the function is zero
249         * @throws MaxIterationsExceededException if the maximum iteration count
250         * is exceeded
251         * @throws FunctionEvaluationException if an error occurs evaluating
252         * the function
253         */
254        private double solve(final UnivariateRealFunction f,
255                             double x0, double y0,
256                             double x1, double y1,
257                             double x2, double y2)
258        throws MaxIterationsExceededException, FunctionEvaluationException {
259    
260            double delta = x1 - x0;
261            double oldDelta = delta;
262    
263            int i = 0;
264            while (i < maximalIterationCount) {
265                if (Math.abs(y2) < Math.abs(y1)) {
266                    // use the bracket point if is better than last approximation
267                    x0 = x1;
268                    x1 = x2;
269                    x2 = x0;
270                    y0 = y1;
271                    y1 = y2;
272                    y2 = y0;
273                }
274                if (Math.abs(y1) <= functionValueAccuracy) {
275                    // Avoid division by very small values. Assume
276                    // the iteration has converged (the problem may
277                    // still be ill conditioned)
278                    setResult(x1, i);
279                    return result;
280                }
281                double dx = x2 - x1;
282                double tolerance =
283                    Math.max(relativeAccuracy * Math.abs(x1), absoluteAccuracy);
284                if (Math.abs(dx) <= tolerance) {
285                    setResult(x1, i);
286                    return result;
287                }
288                if ((Math.abs(oldDelta) < tolerance) ||
289                        (Math.abs(y0) <= Math.abs(y1))) {
290                    // Force bisection.
291                    delta = 0.5 * dx;
292                    oldDelta = delta;
293                } else {
294                    double r3 = y1 / y0;
295                    double p;
296                    double p1;
297                    // the equality test (x0 == x2) is intentional,
298                    // it is part of the original Brent's method,
299                    // it should NOT be replaced by proximity test
300                    if (x0 == x2) {
301                        // Linear interpolation.
302                        p = dx * r3;
303                        p1 = 1.0 - r3;
304                    } else {
305                        // Inverse quadratic interpolation.
306                        double r1 = y0 / y2;
307                        double r2 = y1 / y2;
308                        p = r3 * (dx * r1 * (r1 - r2) - (x1 - x0) * (r2 - 1.0));
309                        p1 = (r1 - 1.0) * (r2 - 1.0) * (r3 - 1.0);
310                    }
311                    if (p > 0.0) {
312                        p1 = -p1;
313                    } else {
314                        p = -p;
315                    }
316                    if (2.0 * p >= 1.5 * dx * p1 - Math.abs(tolerance * p1) ||
317                            p >= Math.abs(0.5 * oldDelta * p1)) {
318                        // Inverse quadratic interpolation gives a value
319                        // in the wrong direction, or progress is slow.
320                        // Fall back to bisection.
321                        delta = 0.5 * dx;
322                        oldDelta = delta;
323                    } else {
324                        oldDelta = delta;
325                        delta = p / p1;
326                    }
327                }
328                // Save old X1, Y1
329                x0 = x1;
330                y0 = y1;
331                // Compute new X1, Y1
332                if (Math.abs(delta) > tolerance) {
333                    x1 = x1 + delta;
334                } else if (dx > 0.0) {
335                    x1 = x1 + 0.5 * tolerance;
336                } else if (dx <= 0.0) {
337                    x1 = x1 - 0.5 * tolerance;
338                }
339                y1 = f.value(x1);
340                if ((y1 > 0) == (y2 > 0)) {
341                    x2 = x0;
342                    y2 = y0;
343                    delta = x1 - x0;
344                    oldDelta = delta;
345                }
346                i++;
347            }
348            throw new MaxIterationsExceededException(maximalIterationCount);
349        }
350    }