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 is based on the false position algorithm.
29 * This implementation is based on Numerical Recipes 3rd ed. Section 9.2
30 * page 451.
31 */
32 public class FalsePositionSingleRootEstimator extends BracketedSingleRootEstimator {
33
34 /**
35 * Maximum allowed number of iterations.
36 */
37 public static final int MAXIT = 30;
38
39 /**
40 * Constant defining default tolerance to find a root.
41 */
42 public static final double DEFAULT_TOLERANCE = 1e-6;
43
44 /**
45 * Minimum allowed tolerance that can be set.
46 */
47 public static final double MIN_TOLERANCE = 0.0;
48
49 /**
50 * Tolerance to find a root. Whenever the variation of the estimated root is
51 * smaller than the provided tolerance, then the algorithm is assumed to be
52 * converged.
53 */
54 private double tolerance;
55
56 /**
57 * Empty constructor.
58 */
59 public FalsePositionSingleRootEstimator() {
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 find a root. During the estimation of a
74 * root, if the variation of the estimated root is smaller than the provided
75 * tolerance, then the algorithm is assumed to be converged, and the root
76 * is ensured to have an accuracy that equals tolerance.
77 * @throws InvalidBracketRangeException Raised if minEvalPoint <
78 * maxEvalPoint.
79 * @throws IllegalArgumentException Raised if provided tolerance is negative.
80 */
81 public FalsePositionSingleRootEstimator(
82 final SingleDimensionFunctionEvaluatorListener listener, final double minEvalPoint,
83 final double maxEvalPoint, final double tolerance) throws InvalidBracketRangeException {
84 super(listener, minEvalPoint, maxEvalPoint);
85 internalSetTolerance(tolerance);
86 }
87
88 /**
89 * Returns tolerance to find a root. Whenever the variation of the estimated
90 * root is smaller than returned tolerance, then the algorithm is assumed to
91 * be converged, and the estimated root is ensured to have an accuracy that
92 * equals the returned tolerance.
93 *
94 * @return Tolerance to find a root.
95 */
96 public double getTolerance() {
97 return tolerance;
98 }
99
100 /**
101 * Sets tolerance to find a root. Whenever the variation of the estimated
102 * root is smaller than provided tolerance, then the algorithm is assumed to
103 * be converged, and the estimated root is ensured to have an accuracy that
104 * equals provided tolerance.
105 *
106 * @param tolerance Tolerance to find a root.
107 * @throws LockedException Raised if this instance is locked while doing
108 * some computations.
109 * @throws IllegalArgumentException Raised if provided tolerance is negative.
110 */
111 public void setTolerance(final double tolerance) throws LockedException {
112 if (isLocked()) {
113 throw new LockedException();
114 }
115 internalSetTolerance(tolerance);
116 }
117
118 /**
119 * Estimates a single root of the provided single dimension function
120 * contained within a given bracket of values.
121 *
122 * @throws LockedException Exception raised if this instance is already
123 * locked.
124 * @throws NotReadyException Exception raised if not enough parameters have
125 * been provided in order to start the estimation.
126 * @throws RootEstimationException Raised if the root estimation failed for
127 * some other reason (usually inability to evaluate the function,
128 * numerical instability or convergence problems, or no roots are found).
129 */
130 @Override
131 public void estimate() throws LockedException, NotReadyException, RootEstimationException {
132 if (isLocked()) {
133 throw new LockedException();
134 }
135 if (!isReady()) {
136 throw new NotReadyException();
137 }
138
139 locked = true;
140 rootAvailable = false;
141
142 try {
143 final var v1 = new double[1];
144 final var v2 = new double[1];
145
146 final var x1 = minEvalPoint;
147 final var x2 = maxEvalPoint;
148 final var xacc = tolerance;
149 double xl;
150 double xh;
151 double del;
152 var fl = listener.evaluate(x1);
153 var fh = listener.evaluate(x2);
154 if (fl * fh > 0.0) {
155 // root must be bracketed
156 locked = false;
157 throw new RootEstimationException();
158 }
159 if (fl < 0.0) {
160 xl = x1;
161 xh = x2;
162 } else {
163 xl = x2;
164 xh = x1;
165
166 v1[0] = fl;
167 v2[0] = fh;
168 swap(v1, v2);
169 fl = v1[0];
170 fh = v2[0];
171 }
172 var dx = xh - xl;
173 for (var j = 0; j < MAXIT; j++) {
174 final var rtf = xl + dx * fl / (fl - fh);
175 final var f = listener.evaluate(rtf);
176 if (f < 0.0) {
177 del = xl - rtf;
178 xl = rtf;
179 fl = f;
180 } else {
181 del = xh - rtf;
182 xh = rtf;
183 fh = f;
184 }
185 dx = xh - xl;
186 if (Math.abs(del) < xacc || f == 0.0) {
187 // result obtained
188 root = rtf;
189 rootAvailable = true;
190 locked = false;
191 return;
192 }
193 }
194 } catch (final EvaluationException e) {
195 throw new RootEstimationException(e);
196 }
197 // too many iterations and error exceeds desired tolerance
198 locked = false;
199 throw new RootEstimationException();
200 }
201
202 /**
203 * Returns boolean indicating whether this instance is ready to compute a
204 * root.
205 * This instance is ready as soon as a listener is provided and a bracket
206 * is provided or computed.
207 *
208 * @return True if instance is ready, false otherwise
209 */
210 @Override
211 public boolean isReady() {
212 return isListenerAvailable() && isBracketAvailable();
213 }
214
215 /**
216 * Internal method to set tolerance to find a root. This method does not
217 * check whether this instance is locked.
218 * Whenever the variation of the estimated root is smaller than provided
219 * tolerance, then the algorithm is assumed to be converged, and the
220 * estimated root is ensured to have an accuracy that equals provided
221 * tolerance.
222 *
223 * @param tolerance Tolerance to find a root.
224 * @throws IllegalArgumentException Raised if provided tolerance is negative.
225 */
226 private void internalSetTolerance(final double tolerance) {
227 if (tolerance < MIN_TOLERANCE) {
228 throw new IllegalArgumentException();
229 }
230 this.tolerance = tolerance;
231 }
232 }