1 /*
2 * Licensed to the Apache Software Foundation (ASF) under one or more
3 * contributor license agreements. See the NOTICE file distributed with
4 * this work for additional information regarding copyright ownership.
5 * The ASF licenses this file to You under the Apache License, Version 2.0
6 * (the "License"); you may not use this file except in compliance with
7 * the License. You may obtain a copy of the License at
8 *
9 * http://www.apache.org/licenses/LICENSE-2.0
10 *
11 * Unless required by applicable law or agreed to in writing, software
12 * distributed under the License is distributed on an "AS IS" BASIS,
13 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14 * See the License for the specific language governing permissions and
15 * limitations under the License.
16 */
17 package org.apache.commons.math3.analysis.solvers;
18
19
20 import org.apache.commons.math3.exception.NoBracketingException;
21 import org.apache.commons.math3.exception.TooManyEvaluationsException;
22 import org.apache.commons.math3.exception.NumberIsTooLargeException;
23 import org.apache.commons.math3.util.FastMath;
24 import org.apache.commons.math3.util.Precision;
25
26 /**
27 * This class implements the <a href="http://mathworld.wolfram.com/BrentsMethod.html">
28 * Brent algorithm</a> for finding zeros of real univariate functions.
29 * The function should be continuous but not necessarily smooth.
30 * The {@code solve} method returns a zero {@code x} of the function {@code f}
31 * in the given interval {@code [a, b]} to within a tolerance
32 * {@code 6 eps abs(x) + t} where {@code eps} is the relative accuracy and
33 * {@code t} is the absolute accuracy.
34 * The given interval must bracket the root.
35 *
36 * @version $Id: BrentSolver.java 1379560 2012-08-31 19:40:30Z erans $
37 */
38 public class BrentSolver extends AbstractUnivariateSolver {
39
40 /** Default absolute accuracy. */
41 private static final double DEFAULT_ABSOLUTE_ACCURACY = 1e-6;
42
43 /**
44 * Construct a solver with default accuracy (1e-6).
45 */
46 public BrentSolver() {
47 this(DEFAULT_ABSOLUTE_ACCURACY);
48 }
49 /**
50 * Construct a solver.
51 *
52 * @param absoluteAccuracy Absolute accuracy.
53 */
54 public BrentSolver(double absoluteAccuracy) {
55 super(absoluteAccuracy);
56 }
57 /**
58 * Construct a solver.
59 *
60 * @param relativeAccuracy Relative accuracy.
61 * @param absoluteAccuracy Absolute accuracy.
62 */
63 public BrentSolver(double relativeAccuracy,
64 double absoluteAccuracy) {
65 super(relativeAccuracy, absoluteAccuracy);
66 }
67 /**
68 * Construct a solver.
69 *
70 * @param relativeAccuracy Relative accuracy.
71 * @param absoluteAccuracy Absolute accuracy.
72 * @param functionValueAccuracy Function value accuracy.
73 */
74 public BrentSolver(double relativeAccuracy,
75 double absoluteAccuracy,
76 double functionValueAccuracy) {
77 super(relativeAccuracy, absoluteAccuracy, functionValueAccuracy);
78 }
79
80 /**
81 * {@inheritDoc}
82 */
83 @Override
84 protected double doSolve()
85 throws NoBracketingException,
86 TooManyEvaluationsException,
87 NumberIsTooLargeException {
88 double min = getMin();
89 double max = getMax();
90 final double initial = getStartValue();
91 final double functionValueAccuracy = getFunctionValueAccuracy();
92
93 verifySequence(min, initial, max);
94
95 // Return the initial guess if it is good enough.
96 double yInitial = computeObjectiveValue(initial);
97 if (FastMath.abs(yInitial) <= functionValueAccuracy) {
98 return initial;
99 }
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 }