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, page
30 * 449
31 */
32 public class SecantSingleRootEstimator extends BracketedSingleRootEstimator {
33
34 /**
35 * Maximum number of iterations.
36 */
37 public static final int MAXIT = 30;
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 SecantSingleRootEstimator() {
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 SecantSingleRootEstimator(
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 * Sets 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 *
103 * @param tolerance Tolerance value.
104 * @throws LockedException Raised if this instance is locked.
105 * @throws IllegalArgumentException Raised if provided tolerance value is
106 * negative.
107 */
108 public void setTolerance(final double tolerance) throws LockedException {
109 if (isLocked()) {
110 throw new LockedException();
111 }
112 internalSetTolerance(tolerance);
113 }
114
115 /**
116 * Estimates a local root for a given single dimension function being
117 * evaluated by provided listener.
118 *
119 * @throws LockedException Exception raised if this instance is already
120 * locked.
121 * @throws NotReadyException Exception raised if either a listener has not
122 * yet been provided or a bracket has not been provided or computed.
123 * @throws RootEstimationException Raised if the root estimation failed for
124 * some other reason (usually inability to evaluate the function,
125 * numerical instability or convergence problems, or no roots are found).
126 */
127 @Override
128 public void estimate() throws LockedException, NotReadyException, RootEstimationException {
129 if (isLocked()) {
130 throw new LockedException();
131 }
132 if (!isReady()) {
133 throw new NotReadyException();
134 }
135
136 locked = true;
137 rootAvailable = false;
138
139 try {
140 final var x1 = minEvalPoint;
141 final var x2 = maxEvalPoint;
142 final var xacc = tolerance;
143 double xl;
144 double rts;
145 var fl = listener.evaluate(x1);
146 var f = listener.evaluate(x2);
147 final var v1 = new double[1];
148 final var v2 = new double[1];
149
150 if (Math.abs(fl) < Math.abs(f)) {
151 rts = x1;
152 xl = x2;
153 v1[0] = fl;
154 v2[0] = f;
155 swap(v1, v2);
156 fl = v1[0];
157 f = v2[0];
158 } else {
159 xl = x1;
160 rts = x2;
161 }
162 for (int j = 0; j < MAXIT; j++) {
163 final var dx = (xl - rts) * f / (f - fl);
164 xl = rts;
165 fl = f;
166 rts += dx;
167 f = listener.evaluate(rts);
168 if (Math.abs(dx) < xacc || f == 0.0) {
169 // Result obtained
170 root = rts;
171 rootAvailable = true;
172 locked = false;
173 return;
174 }
175 }
176 } catch (final EvaluationException e) {
177 throw new RootEstimationException(e);
178 } finally {
179 locked = false;
180 }
181 // too many iterations and error exceeds desired tolerance
182 throw new RootEstimationException();
183 }
184
185 /**
186 * Returns boolean indicating whether this instance is ready to start
187 * estimating a root.
188 * This class will be ready once a listener is provided and a bracket is
189 * either provided or computed.
190 *
191 * @return True if this instance is ready, false otherwise.
192 */
193 @Override
194 public boolean isReady() {
195 return isListenerAvailable() && isBracketAvailable();
196 }
197
198 /**
199 * Internal method to set tolerance value.
200 * Tolerance is the accuracy to be achieved when estimating a root.
201 * If a root is found by this class, it is ensured to have an accuracy below
202 * provided tolerance value.
203 * This method does not check whether this instance is locked or not.
204 *
205 * @param tolerance Tolerance value.
206 * @throws IllegalArgumentException Raised if provided tolerance value is
207 * negative.
208 */
209 private void internalSetTolerance(final double tolerance) {
210 if (tolerance < MIN_TOLERANCE) {
211 throw new IllegalArgumentException();
212 }
213 this.tolerance = tolerance;
214 }
215 }