001 /*
002 * Licensed to the Apache Software Foundation (ASF) under one or more
003 * contributor license agreements. See the NOTICE file distributed with
004 * this work for additional information regarding copyright ownership.
005 * The ASF licenses this file to You under the Apache License, Version 2.0
006 * (the "License"); you may not use this file except in compliance with
007 * the License. You may obtain a copy of the License at
008 *
009 * http://www.apache.org/licenses/LICENSE-2.0
010 *
011 * Unless required by applicable law or agreed to in writing, software
012 * distributed under the License is distributed on an "AS IS" BASIS,
013 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
014 * See the License for the specific language governing permissions and
015 * limitations under the License.
016 */
017 package org.apache.commons.math.analysis.solvers;
018
019
020 import org.apache.commons.math.FunctionEvaluationException;
021 import org.apache.commons.math.MathRuntimeException;
022 import org.apache.commons.math.MaxIterationsExceededException;
023 import org.apache.commons.math.analysis.UnivariateRealFunction;
024
025 /**
026 * Implements the <a href="http://mathworld.wolfram.com/BrentsMethod.html">
027 * Brent algorithm</a> for finding zeros of real univariate functions.
028 * <p>
029 * The function should be continuous but not necessarily smooth.</p>
030 *
031 * @version $Revision:670469 $ $Date:2008-06-23 10:01:38 +0200 (lun., 23 juin 2008) $
032 */
033 public class BrentSolver extends UnivariateRealSolverImpl {
034
035 /**
036 * Default absolute accuracy
037 * @since 2.1
038 */
039 public static final double DEFAULT_ABSOLUTE_ACCURACY = 1E-6;
040
041 /** Default maximum number of iterations
042 * @since 2.1
043 */
044 public static final int DEFAULT_MAXIMUM_ITERATIONS = 100;
045
046 /** Error message for non-bracketing interval. */
047 private static final String NON_BRACKETING_MESSAGE =
048 "function values at endpoints do not have different signs. " +
049 "Endpoints: [{0}, {1}], Values: [{2}, {3}]";
050
051 /** Serializable version identifier */
052 private static final long serialVersionUID = 7694577816772532779L;
053
054 /**
055 * Construct a solver for the given function.
056 *
057 * @param f function to solve.
058 * @deprecated as of 2.0 the function to solve is passed as an argument
059 * to the {@link #solve(UnivariateRealFunction, double, double)} or
060 * {@link UnivariateRealSolverImpl#solve(UnivariateRealFunction, double, double, double)}
061 * method.
062 */
063 @Deprecated
064 public BrentSolver(UnivariateRealFunction f) {
065 super(f, DEFAULT_MAXIMUM_ITERATIONS, DEFAULT_ABSOLUTE_ACCURACY);
066 }
067
068 /**
069 * Construct a solver with default properties.
070 */
071 public BrentSolver() {
072 super(DEFAULT_MAXIMUM_ITERATIONS, DEFAULT_ABSOLUTE_ACCURACY);
073 }
074
075 /**
076 * Construct a solver with the given absolute accuracy.
077 *
078 * @param absoluteAccuracy lower bound for absolute accuracy of solutions returned by the solver
079 * @since 2.1
080 */
081 public BrentSolver(double absoluteAccuracy) {
082 super(DEFAULT_MAXIMUM_ITERATIONS, absoluteAccuracy);
083 }
084
085 /**
086 * Contstruct a solver with the given maximum iterations and absolute accuracy.
087 *
088 * @param maximumIterations maximum number of iterations
089 * @param absoluteAccuracy lower bound for absolute accuracy of solutions returned by the solver
090 * @since 2.1
091 */
092 public BrentSolver(int maximumIterations, double absoluteAccuracy) {
093 super(maximumIterations, absoluteAccuracy);
094 }
095
096 /** {@inheritDoc} */
097 @Deprecated
098 public double solve(double min, double max)
099 throws MaxIterationsExceededException, FunctionEvaluationException {
100 return solve(f, min, max);
101 }
102
103 /** {@inheritDoc} */
104 @Deprecated
105 public double solve(double min, double max, double initial)
106 throws MaxIterationsExceededException, FunctionEvaluationException {
107 return solve(f, min, max, initial);
108 }
109
110 /**
111 * Find a zero in the given interval with an initial guess.
112 * <p>Throws <code>IllegalArgumentException</code> if the values of the
113 * function at the three points have the same sign (note that it is
114 * allowed to have endpoints with the same sign if the initial point has
115 * opposite sign function-wise).</p>
116 *
117 * @param f function to solve.
118 * @param min the lower bound for the interval.
119 * @param max the upper bound for the interval.
120 * @param initial the start value to use (must be set to min if no
121 * initial point is known).
122 * @return the value where the function is zero
123 * @throws MaxIterationsExceededException the maximum iteration count
124 * is exceeded
125 * @throws FunctionEvaluationException if an error occurs evaluating
126 * the function
127 * @throws IllegalArgumentException if initial is not between min and max
128 * (even if it <em>is</em> a root)
129 */
130 public double solve(final UnivariateRealFunction f,
131 final double min, final double max, final double initial)
132 throws MaxIterationsExceededException, FunctionEvaluationException {
133
134 clearResult();
135 if ((initial < min) || (initial > max)) {
136 throw MathRuntimeException.createIllegalArgumentException(
137 "invalid interval, initial value parameters: lower={0}, initial={1}, upper={2}",
138 min, initial, max);
139 }
140
141 // return the initial guess if it is good enough
142 double yInitial = f.value(initial);
143 if (Math.abs(yInitial) <= functionValueAccuracy) {
144 setResult(initial, 0);
145 return result;
146 }
147
148 // return the first endpoint if it is good enough
149 double yMin = f.value(min);
150 if (Math.abs(yMin) <= functionValueAccuracy) {
151 setResult(min, 0);
152 return result;
153 }
154
155 // reduce interval if min and initial bracket the root
156 if (yInitial * yMin < 0) {
157 return solve(f, min, yMin, initial, yInitial, min, yMin);
158 }
159
160 // return the second endpoint if it is good enough
161 double yMax = f.value(max);
162 if (Math.abs(yMax) <= functionValueAccuracy) {
163 setResult(max, 0);
164 return result;
165 }
166
167 // reduce interval if initial and max bracket the root
168 if (yInitial * yMax < 0) {
169 return solve(f, initial, yInitial, max, yMax, initial, yInitial);
170 }
171
172 throw MathRuntimeException.createIllegalArgumentException(
173 NON_BRACKETING_MESSAGE, min, max, yMin, yMax);
174
175 }
176
177 /**
178 * Find a zero in the given interval.
179 * <p>
180 * Requires that the values of the function at the endpoints have opposite
181 * signs. An <code>IllegalArgumentException</code> is thrown if this is not
182 * the case.</p>
183 *
184 * @param f the function to solve
185 * @param min the lower bound for the interval.
186 * @param max the upper bound for the interval.
187 * @return the value where the function is zero
188 * @throws MaxIterationsExceededException if the maximum iteration count is exceeded
189 * @throws FunctionEvaluationException if an error occurs evaluating the
190 * function
191 * @throws IllegalArgumentException if min is not less than max or the
192 * signs of the values of the function at the endpoints are not opposites
193 */
194 public double solve(final UnivariateRealFunction f,
195 final double min, final double max)
196 throws MaxIterationsExceededException,
197 FunctionEvaluationException {
198
199 clearResult();
200 verifyInterval(min, max);
201
202 double ret = Double.NaN;
203
204 double yMin = f.value(min);
205 double yMax = f.value(max);
206
207 // Verify bracketing
208 double sign = yMin * yMax;
209 if (sign > 0) {
210 // check if either value is close to a zero
211 if (Math.abs(yMin) <= functionValueAccuracy) {
212 setResult(min, 0);
213 ret = min;
214 } else if (Math.abs(yMax) <= functionValueAccuracy) {
215 setResult(max, 0);
216 ret = max;
217 } else {
218 // neither value is close to zero and min and max do not bracket root.
219 throw MathRuntimeException.createIllegalArgumentException(
220 NON_BRACKETING_MESSAGE, min, max, yMin, yMax);
221 }
222 } else if (sign < 0){
223 // solve using only the first endpoint as initial guess
224 ret = solve(f, min, yMin, max, yMax, min, yMin);
225 } else {
226 // either min or max is a root
227 if (yMin == 0.0) {
228 ret = min;
229 } else {
230 ret = max;
231 }
232 }
233
234 return ret;
235 }
236
237 /**
238 * Find a zero starting search according to the three provided points.
239 * @param f the function to solve
240 * @param x0 old approximation for the root
241 * @param y0 function value at the approximation for the root
242 * @param x1 last calculated approximation for the root
243 * @param y1 function value at the last calculated approximation
244 * for the root
245 * @param x2 bracket point (must be set to x0 if no bracket point is
246 * known, this will force starting with linear interpolation)
247 * @param y2 function value at the bracket point.
248 * @return the value where the function is zero
249 * @throws MaxIterationsExceededException if the maximum iteration count
250 * is exceeded
251 * @throws FunctionEvaluationException if an error occurs evaluating
252 * the function
253 */
254 private double solve(final UnivariateRealFunction f,
255 double x0, double y0,
256 double x1, double y1,
257 double x2, double y2)
258 throws MaxIterationsExceededException, FunctionEvaluationException {
259
260 double delta = x1 - x0;
261 double oldDelta = delta;
262
263 int i = 0;
264 while (i < maximalIterationCount) {
265 if (Math.abs(y2) < Math.abs(y1)) {
266 // use the bracket point if is better than last approximation
267 x0 = x1;
268 x1 = x2;
269 x2 = x0;
270 y0 = y1;
271 y1 = y2;
272 y2 = y0;
273 }
274 if (Math.abs(y1) <= functionValueAccuracy) {
275 // Avoid division by very small values. Assume
276 // the iteration has converged (the problem may
277 // still be ill conditioned)
278 setResult(x1, i);
279 return result;
280 }
281 double dx = x2 - x1;
282 double tolerance =
283 Math.max(relativeAccuracy * Math.abs(x1), absoluteAccuracy);
284 if (Math.abs(dx) <= tolerance) {
285 setResult(x1, i);
286 return result;
287 }
288 if ((Math.abs(oldDelta) < tolerance) ||
289 (Math.abs(y0) <= Math.abs(y1))) {
290 // Force bisection.
291 delta = 0.5 * dx;
292 oldDelta = delta;
293 } else {
294 double r3 = y1 / y0;
295 double p;
296 double p1;
297 // the equality test (x0 == x2) is intentional,
298 // it is part of the original Brent's method,
299 // it should NOT be replaced by proximity test
300 if (x0 == x2) {
301 // Linear interpolation.
302 p = dx * r3;
303 p1 = 1.0 - r3;
304 } else {
305 // Inverse quadratic interpolation.
306 double r1 = y0 / y2;
307 double r2 = y1 / y2;
308 p = r3 * (dx * r1 * (r1 - r2) - (x1 - x0) * (r2 - 1.0));
309 p1 = (r1 - 1.0) * (r2 - 1.0) * (r3 - 1.0);
310 }
311 if (p > 0.0) {
312 p1 = -p1;
313 } else {
314 p = -p;
315 }
316 if (2.0 * p >= 1.5 * dx * p1 - Math.abs(tolerance * p1) ||
317 p >= Math.abs(0.5 * oldDelta * p1)) {
318 // Inverse quadratic interpolation gives a value
319 // in the wrong direction, or progress is slow.
320 // Fall back to bisection.
321 delta = 0.5 * dx;
322 oldDelta = delta;
323 } else {
324 oldDelta = delta;
325 delta = p / p1;
326 }
327 }
328 // Save old X1, Y1
329 x0 = x1;
330 y0 = y1;
331 // Compute new X1, Y1
332 if (Math.abs(delta) > tolerance) {
333 x1 = x1 + delta;
334 } else if (dx > 0.0) {
335 x1 = x1 + 0.5 * tolerance;
336 } else if (dx <= 0.0) {
337 x1 = x1 - 0.5 * tolerance;
338 }
339 y1 = f.value(x1);
340 if ((y1 > 0) == (y2 > 0)) {
341 x2 = x0;
342 y2 = y0;
343 delta = x1 - x0;
344 oldDelta = delta;
345 }
346 i++;
347 }
348 throw new MaxIterationsExceededException(maximalIterationCount);
349 }
350 }