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 * Finds a single dimensional function's root within a bracket of values using
26 * Newton-Raphson's method.
27 * This class is based on the implementation found in Numerical Recipes 3rd ed.
28 * Section 9.4. page 456.
29 */
30 public class NewtonRaphsonSingleRootEstimator extends DerivativeSingleRootEstimator {
31
32 /**
33 * Maximum number of iterations.
34 */
35 public static final int JMAX = 20;
36
37 /**
38 * Constant defining default accuracy of the estimated root.
39 */
40 public static final double DEFAULT_TOLERANCE = 1e-6;
41
42 /**
43 * Constant defining minimum allowed tolerance.
44 */
45 public static final double MIN_TOLERANCE = 0.0;
46
47 /**
48 * Tolerance value. The algorithm will iterate until the result converges
49 * below this value of accuracy or until the maximum number of iterations is
50 * achieved (and in such case, convergence will be assumed to have failed).
51 */
52 private double tolerance;
53
54 /**
55 * Empty constructor.
56 */
57 public NewtonRaphsonSingleRootEstimator() {
58 super();
59 tolerance = DEFAULT_TOLERANCE;
60 }
61
62 /**
63 * Constructor.
64 *
65 * @param listener Listener to evaluate a single dimension function f(x)
66 * to find its roots.
67 * @param minEvalPoint Smallest value inside the bracket of values where the
68 * root will be searched.
69 * @param maxEvalPoint Largest value inside the bracket of values where the
70 * root will be searched.
71 * @param tolerance Tolerance to be achieved in the estimated root.
72 * @throws InvalidBracketRangeException Raised if minEvalPoint <
73 * maxEvalPoint.
74 * @throws IllegalArgumentException Raised if tolerance is negative.
75 */
76 public NewtonRaphsonSingleRootEstimator(
77 final SingleDimensionFunctionEvaluatorListener listener, final double minEvalPoint,
78 final double maxEvalPoint, final double tolerance) throws InvalidBracketRangeException {
79 super(listener, minEvalPoint, maxEvalPoint);
80 internalSetTolerance(tolerance);
81 }
82
83 /**
84 * Constructor.
85 *
86 * @param listener Listener to evaluate a single dimension function f(x)
87 * to find its roots.
88 * @param derivativeListener Listener to evaluate the function's derivative.
89 * @param minEvalPoint Smallest value inside the bracket of values where the
90 * root will be searched.
91 * @param maxEvalPoint Largest value inside the bracket of values where the
92 * root will be searched.
93 * @param tolerance Tolerance to be achieved in the estimated root.
94 * @throws InvalidBracketRangeException Raised if minEvalPoint <
95 * maxEvalPoint.
96 * @throws IllegalArgumentException Raised if tolerance is negative.
97 */
98 public NewtonRaphsonSingleRootEstimator(
99 final SingleDimensionFunctionEvaluatorListener listener,
100 final SingleDimensionFunctionEvaluatorListener derivativeListener, final double minEvalPoint,
101 final double maxEvalPoint, final double tolerance) throws InvalidBracketRangeException {
102 super(listener, derivativeListener, minEvalPoint, maxEvalPoint);
103 internalSetTolerance(tolerance);
104 }
105
106 /**
107 * Returns tolerance value.
108 * Tolerance is the accuracy to be achieved when estimating a root.
109 * If a root is found by this class, it is ensured to have an accuracy below
110 * the tolerance value.
111 *
112 * @return Tolerance value.
113 */
114 public double getTolerance() {
115 return tolerance;
116 }
117
118 /**
119 * Sets tolerance value.
120 * Tolerance is the accuracy to be achieved when estimating a root.
121 * If a root is found by this class, it is ensured to have an accuracy below
122 * provided tolerance value.
123 *
124 * @param tolerance Tolerance value.
125 * @throws LockedException Raised if this instance is locked.
126 * @throws IllegalArgumentException Raised if provided tolerance value is
127 * negative.
128 */
129 public void setTolerance(final double tolerance) throws LockedException {
130 if (isLocked()) {
131 throw new LockedException();
132 }
133 internalSetTolerance(tolerance);
134 }
135
136 /**
137 * Estimates a local root for a given single dimension function being
138 * evaluated by provided listener.
139 *
140 * @throws LockedException Exception raised if this instance is already
141 * locked.
142 * @throws NotReadyException Exception raised if either a listener has not
143 * yet been provided or a bracket has not been provided or computed.
144 * @throws RootEstimationException Raised if the root estimation failed for
145 * some other reason (usually inability to evaluate the function,
146 * numerical instability or convergence problems, or no roots are found).
147 */
148 @Override
149 @SuppressWarnings("Duplicates")
150 public void estimate() throws LockedException, NotReadyException, RootEstimationException {
151
152 if (isLocked()) {
153 throw new LockedException();
154 }
155 if (!isReady()) {
156 throw new NotReadyException();
157 }
158
159 locked = true;
160 rootAvailable = false;
161
162 final var x1 = minEvalPoint;
163 final var x2 = maxEvalPoint;
164 final var xacc = tolerance;
165
166 var rtn = 0.5 * (x1 + x2);
167 double f;
168 double df;
169 for (var j = 0; j < JMAX; j++) {
170 try {
171 f = listener.evaluate(rtn);
172 df = derivativeListener.evaluate(rtn);
173 } catch (final EvaluationException e) {
174 throw new RootEstimationException(e);
175 }
176
177 final var dx = f / df;
178 rtn -= dx;
179 if ((x1 - rtn) * (rtn - x2) < 0.0) {
180 // jumped out of brackets
181 locked = false;
182 throw new RootEstimationException();
183 }
184 if (Math.abs(dx) < xacc) {
185 // root found
186 root = rtn;
187 rootAvailable = true;
188 locked = false;
189 return;
190 }
191 }
192 // maximum number of iterations exceeded
193 locked = false;
194 throw new RootEstimationException();
195 }
196
197 /**
198 * Internal method to set tolerance value.
199 * Tolerance is the accuracy to be achieved when estimating a root.
200 * If a root is found by this class, it is ensured to have an accuracy below
201 * provided tolerance value.
202 * This method does not check whether this instance is locked or not.
203 *
204 * @param tolerance Tolerance value.
205 * @throws IllegalArgumentException Raised if provided tolerance value is
206 * negative.
207 */
208 private void internalSetTolerance(final double tolerance) {
209 if (tolerance < MIN_TOLERANCE) {
210 throw new IllegalArgumentException();
211 }
212 this.tolerance = tolerance;
213 }
214 }