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