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.math3.analysis.solvers;
018
019
020 import org.apache.commons.math3.exception.NoBracketingException;
021 import org.apache.commons.math3.exception.TooManyEvaluationsException;
022 import org.apache.commons.math3.exception.NumberIsTooLargeException;
023 import org.apache.commons.math3.util.FastMath;
024 import 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 */
038 public 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 }