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.EvaluationException;
19 import com.irurueta.numerical.InvalidBracketRangeException;
20 import com.irurueta.numerical.LockedException;
21 import com.irurueta.numerical.NotReadyException;
22 import com.irurueta.numerical.SingleDimensionFunctionEvaluatorListener;
23
24 /**
25 * Computes a root for a single dimension function inside a given bracket of
26 * values, in other words, root will only be searched within provided minimum
27 * and maximum evaluation points.
28 * This class searches for REAL roots only!
29 * This implementation is based on Numerical Recipes 3rd ed. Section 9.2.1, page
30 * 452.
31 */
32 public class RidderSingleRootEstimator extends BracketedSingleRootEstimator {
33
34 /**
35 * Maximum number of iterations.
36 */
37 public static final int MAXIT = 60;
38
39 /**
40 * Constant defining default accuracy of the estimated root.
41 */
42 public static final double DEFAULT_TOLERANCE = 1e-6;
43
44 /**
45 * Constant defining minimum allowed tolerance.
46 */
47 public static final double MIN_TOLERANCE = 0.0;
48
49 /**
50 * Tolerance value. The algorithm will iterate until the result converges
51 * below this value of accuracy or until the maximum number of iterations is
52 * achieved (and in such case, convergence will be assumed to have failed).
53 */
54 private double tolerance;
55
56 /**
57 * Empty constructor.
58 */
59 public RidderSingleRootEstimator() {
60 super();
61 tolerance = DEFAULT_TOLERANCE;
62 }
63
64 /**
65 * Constructor.
66 *
67 * @param listener Listener to evaluate a single dimension function f(x)
68 * to find its roots.
69 * @param minEvalPoint Smallest value inside the bracket of values where the
70 * root will be searched.
71 * @param maxEvalPoint Largest value inside the bracket of values where the
72 * root will be searched.
73 * @param tolerance Tolerance to be achieved in the estimated root.
74 * @throws InvalidBracketRangeException Raised if minEvalPoint <
75 * maxEvalPoint.
76 * @throws IllegalArgumentException Raised if tolerance is negative.
77 */
78 public RidderSingleRootEstimator(
79 final SingleDimensionFunctionEvaluatorListener listener, final double minEvalPoint,
80 final double maxEvalPoint, final double tolerance) throws InvalidBracketRangeException {
81 super(listener, minEvalPoint, maxEvalPoint);
82 internalSetTolerance(tolerance);
83 }
84
85 /**
86 * Returns tolerance value.
87 * Tolerance is the accuracy to be achieved when estimating a root.
88 * If a root is found by this class, it is ensured to have an accuracy below
89 * the tolerance value.
90 *
91 * @return Tolerance value.
92 */
93 public double getTolerance() {
94 return tolerance;
95 }
96
97 /**
98 * Internal method to set tolerance value.
99 * Tolerance is the accuracy to be achieved when estimating a root.
100 * If a root is found by this class, it is ensured to have an accuracy below
101 * provided tolerance value.
102 * This method does not check whether this instance is locked or not.
103 *
104 * @param tolerance Tolerance value.
105 * @throws IllegalArgumentException Raised if provided tolerance value is
106 * negative.
107 */
108 private void internalSetTolerance(final double tolerance) {
109 if (tolerance < MIN_TOLERANCE) {
110 throw new IllegalArgumentException();
111 }
112 this.tolerance = tolerance;
113 }
114
115 /**
116 * Sets tolerance value.
117 * Tolerance is the accuracy to be achieved when estimating a root.
118 * If a root is found by this class, it is ensured to have an accuracy below
119 * provided tolerance value.
120 *
121 * @param tolerance Tolerance value.
122 * @throws LockedException Raised if this instance is locked.
123 * @throws IllegalArgumentException Raised if provided tolerance value is
124 * negative.
125 */
126 public void setTolerance(final double tolerance) throws LockedException {
127 if (isLocked()) {
128 throw new LockedException();
129 }
130 internalSetTolerance(tolerance);
131 }
132
133 /**
134 * Estimates a local root for a given single dimension function being
135 * evaluated by provided listener.
136 *
137 * @throws LockedException Exception raised if this instance is already
138 * locked.
139 * @throws NotReadyException Exception raised if either a listener has not
140 * yet been provided or a bracket has not been provided or computed.
141 * @throws RootEstimationException Raised if the root estimation failed for
142 * some other reason (usually inability to evaluate the function,
143 * numerical instability or convergence problems, or no roots are found).
144 */
145 @Override
146 @SuppressWarnings("Duplicates")
147 public void estimate() throws LockedException, NotReadyException, RootEstimationException {
148 if (isLocked()) {
149 throw new LockedException();
150 }
151 if (!isReady()) {
152 throw new NotReadyException();
153 }
154
155 locked = true;
156
157 final var x1 = minEvalPoint;
158 final var x2 = maxEvalPoint;
159 final var xacc = tolerance;
160 double fl;
161 double fh;
162 try {
163 fl = listener.evaluate(x1);
164 fh = listener.evaluate(x2);
165 } catch (final EvaluationException e) {
166 throw new RootEstimationException(e);
167 }
168
169 double ans;
170 var found = false;
171 if ((fl > 0.0 && fh < 0.0) || (fl < 0.0 && fh > 0.0)) {
172 var xl = x1;
173 var xh = x2;
174 ans = -9.99e99;
175 try {
176 for (var j = 0; j < MAXIT; j++) {
177 final var xm = 0.5 * (xl + xh);
178 final var fm = listener.evaluate(xm);
179 final var s = Math.sqrt(fm * fm - fl * fh);
180 if (s == 0.0) {
181 found = true;
182 break;
183 }
184 final var xnew = xm + (xm - xl) * ((fl >= fh ? 1.0 : -1.0) * fm / s);
185 if (Math.abs(xnew - ans) <= xacc) {
186 // result found
187 found = true;
188 break;
189 }
190 ans = xnew;
191 final var fnew = listener.evaluate(ans);
192 if (sign(fm, fnew) != fm) {
193 xl = xm;
194 fl = fm;
195 xh = ans;
196 fh = fnew;
197 } else if (sign(fl, fnew) != fl) {
198 xh = ans;
199 fh = fnew;
200 } else if (sign(fh, fnew) != fh) {
201 xl = ans;
202 fl = fnew;
203 } else {
204 // never get here
205 locked = false;
206 throw new RootEstimationException();
207 }
208 if (Math.abs(xh - xl) <= xacc) {
209 // result found
210 found = true;
211 break;
212 }
213
214 }
215 } catch (final EvaluationException e) {
216 throw new RootEstimationException(e);
217 }
218 if (!found) {
219 // too many iterations and error exceeds desired tolerance
220 locked = false;
221 throw new RootEstimationException();
222 }
223 } else {
224 if (fl == 0.0) {
225 // result found
226 ans = x1;
227 } else if (fh == 0.0) {
228 // result found
229 ans = x2;
230 } else {
231 locked = false;
232 throw new RootEstimationException();
233 }
234 }
235
236 // result found
237 root = ans;
238 rootAvailable = true;
239 locked = false;
240 }
241
242 /**
243 * Returns boolean indicating whether this instance is ready to start
244 * estimating a root.
245 * This class will be ready once a listener is provided and a bracket is
246 * either provided or computed.
247 *
248 * @return True if this instance is ready, false otherwise.
249 */
250 @Override
251 public boolean isReady() {
252 return isListenerAvailable() && isBracketAvailable();
253 }
254 }