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.numbers.rootfinder;
18
19 import java.util.function.DoubleUnaryOperator;
20
21 /**
22 * This class implements the <a href="http://mathworld.wolfram.com/BrentsMethod.html">
23 * Brent algorithm</a> for finding zeros of real univariate functions.
24 * The function should be continuous but not necessarily smooth.
25 * The {@code solve} method returns a zero {@code x} of the function {@code f}
26 * in the given interval {@code [a, b]} to within a tolerance
27 * {@code 2 eps abs(x) + t} where {@code eps} is the relative accuracy and
28 * {@code t} is the absolute accuracy.
29 * <p>The given interval must bracket the root.</p>
30 * <p>
31 * The reference implementation is given in chapter 4 of
32 * <blockquote>
33 * <b>Algorithms for Minimization Without Derivatives</b>,
34 * <em>Richard P. Brent</em>,
35 * Dover, 2002
36 * </blockquote>
37 */
38 public class BrentSolver {
39 /** Relative accuracy. */
40 private final double relativeAccuracy;
41 /** Absolute accuracy. */
42 private final double absoluteAccuracy;
43 /** Function accuracy. */
44 private final double functionValueAccuracy;
45
46 /**
47 * Construct a solver.
48 *
49 * @param relativeAccuracy Relative accuracy.
50 * @param absoluteAccuracy Absolute accuracy.
51 * @param functionValueAccuracy Function value accuracy.
52 */
53 public BrentSolver(double relativeAccuracy,
54 double absoluteAccuracy,
55 double functionValueAccuracy) {
56 this.relativeAccuracy = relativeAccuracy;
57 this.absoluteAccuracy = absoluteAccuracy;
58 this.functionValueAccuracy = functionValueAccuracy;
59 }
60
61 /**
62 * Search the function's zero within the given interval.
63 *
64 * @param func Function to solve.
65 * @param min Lower bound.
66 * @param max Upper bound.
67 * @return the root.
68 * @throws IllegalArgumentException if {@code min > max}.
69 * @throws IllegalArgumentException if the given interval does
70 * not bracket the root.
71 */
72 public double findRoot(DoubleUnaryOperator func,
73 double min,
74 double max) {
75 // Avoid overflow computing the initial value: 0.5 * (min + max)
76 // Note: This sum is invalid if min == max == Double.MIN_VALUE
77 // so detect this edge case. It will raise a bracketing exception
78 // if min is not the root within the configured function accuracy;
79 // otherwise min is returned.
80 final double initial = min == max ? min : 0.5 * min + 0.5 * max;
81 return findRoot(func, min, initial, max);
82 }
83
84 /**
85 * Search the function's zero within the given interval,
86 * starting from the given estimate.
87 *
88 * @param func Function to solve.
89 * @param min Lower bound.
90 * @param initial Initial guess.
91 * @param max Upper bound.
92 * @return the root.
93 * @throws IllegalArgumentException if {@code min > max} or
94 * {@code initial} is not in the {@code [min, max]} interval.
95 * @throws IllegalArgumentException if the given interval does
96 * not bracket the root.
97 */
98 public double findRoot(DoubleUnaryOperator func,
99 double min,
100 double initial,
101 double max) {
102 if (min > max) {
103 throw new SolverException(SolverException.TOO_LARGE, min, max);
104 }
105 if (initial < min ||
106 initial > max) {
107 throw new SolverException(SolverException.OUT_OF_RANGE, initial, min, max);
108 }
109
110 // Return the initial guess if it is good enough.
111 final double yInitial = func.applyAsDouble(initial);
112 if (Math.abs(yInitial) <= functionValueAccuracy) {
113 return initial;
114 }
115
116 // Return the first endpoint if it is good enough.
117 final double yMin = func.applyAsDouble(min);
118 if (Math.abs(yMin) <= functionValueAccuracy) {
119 return min;
120 }
121
122 // Reduce interval if min and initial bracket the root.
123 if (Double.compare(yInitial * yMin, 0.0) < 0) {
124 return brent(func, min, initial, yMin, yInitial);
125 }
126
127 // Return the second endpoint if it is good enough.
128 final double yMax = func.applyAsDouble(max);
129 if (Math.abs(yMax) <= functionValueAccuracy) {
130 return max;
131 }
132
133 // Reduce interval if initial and max bracket the root.
134 if (Double.compare(yInitial * yMax, 0.0) < 0) {
135 return brent(func, initial, max, yInitial, yMax);
136 }
137
138 throw new SolverException(SolverException.BRACKETING, min, yMin, max, yMax);
139 }
140
141 /**
142 * Search for a zero inside the provided interval.
143 * This implementation is based on the algorithm described at page 58 of
144 * the book
145 * <blockquote>
146 * <b>Algorithms for Minimization Without Derivatives</b>,
147 * <i>Richard P. Brent</i>,
148 * Dover 0-486-41998-3
149 * </blockquote>
150 *
151 * @param func Function to solve.
152 * @param lo Lower bound of the search interval.
153 * @param hi Higher bound of the search interval.
154 * @param fLo Function value at the lower bound of the search interval.
155 * @param fHi Function value at the higher bound of the search interval.
156 * @return the value where the function is zero.
157 */
158 private double brent(DoubleUnaryOperator func,
159 double lo, double hi,
160 double fLo, double fHi) {
161 double a = lo;
162 double fa = fLo;
163 double b = hi;
164 double fb = fHi;
165 double c = a;
166 double fc = fa;
167 double d = b - a;
168 double e = d;
169
170 final double t = absoluteAccuracy;
171 final double eps = relativeAccuracy;
172
173 while (true) {
174 if (Math.abs(fc) < Math.abs(fb)) {
175 a = b;
176 b = c;
177 c = a;
178 fa = fb;
179 fb = fc;
180 fc = fa;
181 }
182
183 final double tol = 2 * eps * Math.abs(b) + t;
184 final double m = 0.5 * (c - b);
185
186 if (Math.abs(m) <= tol ||
187 equalsZero(fb)) {
188 return b;
189 }
190 if (Math.abs(e) < tol ||
191 Math.abs(fa) <= Math.abs(fb)) {
192 // Force bisection.
193 d = m;
194 e = d;
195 } else {
196 final double s = fb / fa;
197 double p;
198 double q;
199 // The equality test (a == c) is intentional,
200 // it is part of the original Brent's method and
201 // it should NOT be replaced by proximity test.
202 if (a == c) {
203 // Linear interpolation.
204 p = 2 * m * s;
205 q = 1 - s;
206 } else {
207 // Inverse quadratic interpolation.
208 q = fa / fc;
209 final double r = fb / fc;
210 p = s * (2 * m * q * (q - r) - (b - a) * (r - 1));
211 q = (q - 1) * (r - 1) * (s - 1);
212 }
213 if (p > 0) {
214 q = -q;
215 } else {
216 p = -p;
217 }
218 if (p >= 1.5 * m * q - Math.abs(tol * q) ||
219 p >= Math.abs(0.5 * e * q)) {
220 // Inverse quadratic interpolation gives a value
221 // in the wrong direction, or progress is slow.
222 // Fall back to bisection.
223 d = m;
224 e = d;
225 } else {
226 e = d;
227 d = p / q;
228 }
229 }
230 a = b;
231 fa = fb;
232
233 if (Math.abs(d) > tol) {
234 b += d;
235 } else if (m > 0) {
236 b += tol;
237 } else {
238 b -= tol;
239 }
240 fb = func.applyAsDouble(b);
241 if ((fb > 0 && fc > 0) ||
242 (fb <= 0 && fc <= 0)) {
243 c = a;
244 fc = fa;
245 d = b - a;
246 e = d;
247 }
248 }
249 }
250
251 /**
252 * Return true if the value is within 1 ULP of zero.
253 *
254 * @param value Value
255 * @return true if zero within a 1 ULP tolerance
256 */
257 private static boolean equalsZero(double value) {
258 return Math.abs(value) <= Double.MIN_VALUE;
259 }
260 }