1 /*
2 * Copyright (C) 2012 Alberto Irurueta Carro (alberto@irurueta.com)
3 *
4 * Licensed under the Apache License, Version 2.0 (the "License");
5 * you may not use this file except in compliance with the License.
6 * You may obtain a copy of the License at
7 *
8 * http://www.apache.org/licenses/LICENSE-2.0
9 *
10 * Unless required by applicable law or agreed to in writing, software
11 * distributed under the License is distributed on an "AS IS" BASIS,
12 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
13 * See the License for the specific language governing permissions and
14 * limitations under the License.
15 */
16 package com.irurueta.numerical.roots;
17
18 import com.irurueta.numerical.*;
19
20 /**
21 * This class estimates the root of a single dimension continuous function using
22 * Brent's method.
23 * The implementation of this class is based on Numerical Recipes 3rd ed.
24 * Section 9.3. Page 454.
25 */
26 public class BrentSingleRootEstimator extends BracketedSingleRootEstimator {
27
28 /**
29 * Constant defining maximum number of iterations.
30 */
31 public static final int ITMAX = 100;
32
33 /**
34 * Constant defining a small value which is considered as machine precision.
35 */
36 public static final double EPS = 1e-10;
37
38 /**
39 * Constant defining default accuracy of the estimated root.
40 */
41 public static final double DEFAULT_TOLERANCE = 1e-6;
42
43 /**
44 * Constant defining minimum allowed tolerance.
45 */
46 public static final double MIN_TOLERANCE = 0.0;
47
48 /**
49 * Tolerance value. The algorithm will iterate until the result converges
50 * below this value of accuracy or until the maximum number of iterations is
51 * achieved (and in such case, convergence will be assumed to have failed).
52 */
53 private double tolerance;
54
55 /**
56 * Empty constructor.
57 */
58 public BrentSingleRootEstimator() {
59 super();
60 tolerance = DEFAULT_TOLERANCE;
61 }
62
63 /**
64 * Constructor.
65 *
66 * @param listener Listener to evaluate a single dimension function f(x)
67 * to find its roots.
68 * @param minEvalPoint Smallest value inside the bracket of values where the
69 * root will be searched.
70 * @param maxEvalPoint Largest value inside the bracket of values where the
71 * root will be searched.
72 * @param tolerance Tolerance to be achieved in the estimated root.
73 * @throws InvalidBracketRangeException Raised if minEvalPoint <
74 * maxEvalPoint.
75 * @throws IllegalArgumentException Raised if tolerance is negative.
76 */
77 public BrentSingleRootEstimator(
78 final SingleDimensionFunctionEvaluatorListener listener, final double minEvalPoint,
79 final double maxEvalPoint, final double tolerance) throws InvalidBracketRangeException {
80 super(listener, minEvalPoint, maxEvalPoint);
81 internalSetTolerance(tolerance);
82 }
83
84 /**
85 * Returns tolerance value.
86 * Tolerance is the accuracy to be achieved when estimating a root.
87 * If a root is found by this class, it is ensured to have an accuracy below
88 * the tolerance value.
89 *
90 * @return Tolerance value.
91 */
92 public double getTolerance() {
93 return tolerance;
94 }
95
96 /**
97 * Internal method to set tolerance value.
98 * Tolerance is the accuracy to be achieved when estimating a root.
99 * If a root is found by this class, it is ensured to have an accuracy below
100 * provided tolerance value.
101 * This method does not check whether this instance is locked or not.
102 *
103 * @param tolerance Tolerance value.
104 * @throws IllegalArgumentException Raised if provided tolerance value is
105 * negative.
106 */
107 private void internalSetTolerance(final double tolerance) {
108 if (tolerance < MIN_TOLERANCE) {
109 throw new IllegalArgumentException();
110 }
111 this.tolerance = tolerance;
112 }
113
114 /**
115 * Sets tolerance value.
116 * Tolerance is the accuracy to be achieved when estimating a root.
117 * If a root is found by this class, it is ensured to have an accuracy below
118 * provided tolerance value.
119 *
120 * @param tolerance Tolerance value.
121 * @throws LockedException Raised if this instance is locked.
122 * @throws IllegalArgumentException Raised if provided tolerance value is
123 * negative.
124 */
125 public void setTolerance(final double tolerance) throws LockedException {
126 if (isLocked()) {
127 throw new LockedException();
128 }
129 internalSetTolerance(tolerance);
130 }
131
132 /**
133 * Estimates a local root for a given single dimension function being
134 * evaluated by provided listener.
135 *
136 * @throws LockedException Exception raised if this instance is already
137 * locked.
138 * @throws NotReadyException Exception raised if either a listener has not
139 * yet been provided or a bracket has not been provided or computed.
140 * @throws RootEstimationException Raised if the root estimation failed for
141 * some other reason (usually inability to evaluate the function,
142 * numerical instability or convergence problems, or no roots are found).
143 */
144 @Override
145 @SuppressWarnings("Duplicates")
146 public void estimate() throws LockedException, NotReadyException, RootEstimationException {
147 if (isLocked()) {
148 throw new LockedException();
149 }
150 if (!isReady()) {
151 throw new NotReadyException();
152 }
153
154 locked = true;
155
156 final var x1 = minEvalPoint;
157 final var x2 = maxEvalPoint;
158 final var tol = tolerance;
159 var a = x1;
160 var b = x2;
161 var c = x2;
162 var d = 0.0;
163 var e = 0.0;
164 double fc;
165 double p;
166 double q;
167 double r;
168 double s;
169 double tol1;
170 double xm;
171 double fa;
172 double fb;
173 try {
174 fa = listener.evaluate(a);
175 fb = listener.evaluate(b);
176 } catch (final EvaluationException ex) {
177 throw new RootEstimationException(ex);
178 }
179
180 if ((fa > 0.0 && fb > 0.0) || (fa < 0.0 && fb < 0.0)) {
181 // root must be bracketed
182 locked = false;
183 throw new RootEstimationException();
184 }
185 fc = fb;
186 for (var iter = 0; iter < ITMAX; iter++) {
187 if ((fb > 0.0 && fc > 0.0) || (fb < 0.0 && fc < 0.0)) {
188 c = a;
189 fc = fa;
190 e = d = b - a;
191 }
192 if (Math.abs(fc) < Math.abs(fb)) {
193 a = b;
194 b = c;
195 c = a;
196 fa = fb;
197 fb = fc;
198 fc = fa;
199 }
200 tol1 = 2.0 * EPS * Math.abs(b) + 0.5 * tol;
201 xm = 0.5 * (c - b);
202 if (Math.abs(xm) <= tol1 || fb == 0.0) {
203 // root found
204 root = b;
205 rootAvailable = true;
206 locked = false;
207 return;
208 }
209 if (Math.abs(e) >= tol1 && Math.abs(fa) > Math.abs(fb)) {
210 s = fb / fa;
211 if (a == c) {
212 p = 2.0 * xm * s;
213 q = 1.0 - s;
214 } else {
215 q = fa / fc;
216 r = fb / fc;
217 p = s * (2.0 * xm * q * (q - r) - (b - a) * (r - 1.0));
218 q = (q - 1.0) * (r - 1.0) * (s - 1.0);
219 }
220 if (p > 0.0) {
221 q = -q;
222 }
223 p = Math.abs(p);
224 final double min1 = 3.0 * xm * q - Math.abs(tol1 * q);
225 final double min2 = Math.abs(e * q);
226 if (2.0 * p < (Math.min(min1, min2))) {
227 e = d;
228 d = p / q;
229 } else {
230 d = xm;
231 e = d;
232 }
233 } else {
234 d = xm;
235 e = d;
236 }
237 a = b;
238 fa = fb;
239 if (Math.abs(d) > tol1) {
240 b += d;
241 } else {
242 b += sign(tol1, xm);
243 }
244 try {
245 fb = listener.evaluate(b);
246 } catch (final EvaluationException ex) {
247 throw new RootEstimationException(ex);
248 }
249 }
250 // maximum number of iterations exceeded
251 locked = false;
252 throw new RootEstimationException();
253 }
254
255 /**
256 * Returns boolean indicating whether this instance is ready to start
257 * estimating a root.
258 * This class will be ready once a listener is provided and a bracket is
259 * either provided or computed.
260 *
261 * @return True if this instance is ready, false otherwise.
262 */
263 @Override
264 public boolean isReady() {
265 return isListenerAvailable() && isBracketAvailable();
266 }
267 }