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 */
017package org.apache.commons.math3.analysis.solvers;
018
019
020import org.apache.commons.math3.exception.NoBracketingException;
021import org.apache.commons.math3.exception.TooManyEvaluationsException;
022import org.apache.commons.math3.exception.NumberIsTooLargeException;
023import org.apache.commons.math3.util.FastMath;
024import org.apache.commons.math3.util.Precision;
025
026/**
027 * This class implements the <a href="http://mathworld.wolfram.com/BrentsMethod.html">
028 * Brent algorithm</a> for finding zeros of real univariate functions.
029 * The function should be continuous but not necessarily smooth.
030 * The {@code solve} method returns a zero {@code x} of the function {@code f}
031 * in the given interval {@code [a, b]} to within a tolerance
032 * {@code 2 eps abs(x) + t} where {@code eps} is the relative accuracy and
033 * {@code t} is the absolute accuracy.
034 * The given interval must bracket the root.
035 * <p>
036 *  The reference implementation is given in chapter 4 of
037 *  <blockquote>
038 *   <b>Algorithms for Minimization Without Derivatives</b><br>
039 *   <em>Richard P. Brent</em><br>
040 *   Dover, 2002<br>
041 *  </blockquote>
042 * </p>
043 *
044 * @see BaseAbstractUnivariateSolver
045 */
046public class BrentSolver extends AbstractUnivariateSolver {
047
048    /** Default absolute accuracy. */
049    private static final double DEFAULT_ABSOLUTE_ACCURACY = 1e-6;
050
051    /**
052     * Construct a solver with default absolute accuracy (1e-6).
053     */
054    public BrentSolver() {
055        this(DEFAULT_ABSOLUTE_ACCURACY);
056    }
057    /**
058     * Construct a solver.
059     *
060     * @param absoluteAccuracy Absolute accuracy.
061     */
062    public BrentSolver(double absoluteAccuracy) {
063        super(absoluteAccuracy);
064    }
065    /**
066     * Construct a solver.
067     *
068     * @param relativeAccuracy Relative accuracy.
069     * @param absoluteAccuracy Absolute accuracy.
070     */
071    public BrentSolver(double relativeAccuracy,
072                       double absoluteAccuracy) {
073        super(relativeAccuracy, absoluteAccuracy);
074    }
075    /**
076     * Construct a solver.
077     *
078     * @param relativeAccuracy Relative accuracy.
079     * @param absoluteAccuracy Absolute accuracy.
080     * @param functionValueAccuracy Function value accuracy.
081     *
082     * @see BaseAbstractUnivariateSolver#BaseAbstractUnivariateSolver(double,double,double)
083     */
084    public BrentSolver(double relativeAccuracy,
085                       double absoluteAccuracy,
086                       double functionValueAccuracy) {
087        super(relativeAccuracy, absoluteAccuracy, functionValueAccuracy);
088    }
089
090    /**
091     * {@inheritDoc}
092     */
093    @Override
094    protected double doSolve()
095        throws NoBracketingException,
096               TooManyEvaluationsException,
097               NumberIsTooLargeException {
098        double min = getMin();
099        double max = getMax();
100        final double initial = getStartValue();
101        final double functionValueAccuracy = getFunctionValueAccuracy();
102
103        verifySequence(min, initial, max);
104
105        // Return the initial guess if it is good enough.
106        double yInitial = computeObjectiveValue(initial);
107        if (FastMath.abs(yInitial) <= functionValueAccuracy) {
108            return initial;
109        }
110
111        // Return the first endpoint if it is good enough.
112        double yMin = computeObjectiveValue(min);
113        if (FastMath.abs(yMin) <= functionValueAccuracy) {
114            return min;
115        }
116
117        // Reduce interval if min and initial bracket the root.
118        if (yInitial * yMin < 0) {
119            return brent(min, initial, yMin, yInitial);
120        }
121
122        // Return the second endpoint if it is good enough.
123        double yMax = computeObjectiveValue(max);
124        if (FastMath.abs(yMax) <= functionValueAccuracy) {
125            return max;
126        }
127
128        // Reduce interval if initial and max bracket the root.
129        if (yInitial * yMax < 0) {
130            return brent(initial, max, yInitial, yMax);
131        }
132
133        throw new NoBracketingException(min, max, yMin, yMax);
134    }
135
136    /**
137     * Search for a zero inside the provided interval.
138     * This implementation is based on the algorithm described at page 58 of
139     * the book
140     * <blockquote>
141     *  <b>Algorithms for Minimization Without Derivatives</b>
142     *  <it>Richard P. Brent</it>
143     *  Dover 0-486-41998-3
144     * </blockquote>
145     *
146     * @param lo Lower bound of the search interval.
147     * @param hi Higher bound of the search interval.
148     * @param fLo Function value at the lower bound of the search interval.
149     * @param fHi Function value at the higher bound of the search interval.
150     * @return the value where the function is zero.
151     */
152    private double brent(double lo, double hi,
153                         double fLo, double fHi) {
154        double a = lo;
155        double fa = fLo;
156        double b = hi;
157        double fb = fHi;
158        double c = a;
159        double fc = fa;
160        double d = b - a;
161        double e = d;
162
163        final double t = getAbsoluteAccuracy();
164        final double eps = getRelativeAccuracy();
165
166        while (true) {
167            if (FastMath.abs(fc) < FastMath.abs(fb)) {
168                a = b;
169                b = c;
170                c = a;
171                fa = fb;
172                fb = fc;
173                fc = fa;
174            }
175
176            final double tol = 2 * eps * FastMath.abs(b) + t;
177            final double m = 0.5 * (c - b);
178
179            if (FastMath.abs(m) <= tol ||
180                Precision.equals(fb, 0))  {
181                return b;
182            }
183            if (FastMath.abs(e) < tol ||
184                FastMath.abs(fa) <= FastMath.abs(fb)) {
185                // Force bisection.
186                d = m;
187                e = d;
188            } else {
189                double s = fb / fa;
190                double p;
191                double q;
192                // The equality test (a == c) is intentional,
193                // it is part of the original Brent's method and
194                // it should NOT be replaced by proximity test.
195                if (a == c) {
196                    // Linear interpolation.
197                    p = 2 * m * s;
198                    q = 1 - s;
199                } else {
200                    // Inverse quadratic interpolation.
201                    q = fa / fc;
202                    final double r = fb / fc;
203                    p = s * (2 * m * q * (q - r) - (b - a) * (r - 1));
204                    q = (q - 1) * (r - 1) * (s - 1);
205                }
206                if (p > 0) {
207                    q = -q;
208                } else {
209                    p = -p;
210                }
211                s = e;
212                e = d;
213                if (p >= 1.5 * m * q - FastMath.abs(tol * q) ||
214                    p >= FastMath.abs(0.5 * s * q)) {
215                    // Inverse quadratic interpolation gives a value
216                    // in the wrong direction, or progress is slow.
217                    // Fall back to bisection.
218                    d = m;
219                    e = d;
220                } else {
221                    d = p / q;
222                }
223            }
224            a = b;
225            fa = fb;
226
227            if (FastMath.abs(d) > tol) {
228                b += d;
229            } else if (m > 0) {
230                b += tol;
231            } else {
232                b -= tol;
233            }
234            fb = computeObjectiveValue(b);
235            if ((fb > 0 && fc > 0) ||
236                (fb <= 0 && fc <= 0)) {
237                c = a;
238                fc = fa;
239                d = b - a;
240                e = d;
241            }
242        }
243    }
244}