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