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 /** Serializable version identifier */
036 private static final long serialVersionUID = 7694577816772532779L;
037
038 /**
039 * Construct a solver for the given function.
040 *
041 * @param f function to solve.
042 * @deprecated as of 2.0 the function to solve is passed as an argument
043 * to the {@link #solve(UnivariateRealFunction, double, double)} or
044 * {@link UnivariateRealSolverImpl#solve(UnivariateRealFunction, double, double, double)}
045 * method.
046 */
047 @Deprecated
048 public BrentSolver(UnivariateRealFunction f) {
049 super(f, 100, 1E-6);
050 }
051
052 /**
053 * Construct a solver.
054 */
055 public BrentSolver() {
056 super(100, 1E-6);
057 }
058
059 /** {@inheritDoc} */
060 @Deprecated
061 public double solve(double min, double max)
062 throws MaxIterationsExceededException, FunctionEvaluationException {
063 return solve(f, min, max);
064 }
065
066 /** {@inheritDoc} */
067 @Deprecated
068 public double solve(double min, double max, double initial)
069 throws MaxIterationsExceededException, FunctionEvaluationException {
070 return solve(f, min, max, initial);
071 }
072
073 /**
074 * Find a zero in the given interval with an initial guess.
075 * <p>Throws <code>IllegalArgumentException</code> if the values of the
076 * function at the three points have the same sign (note that it is
077 * allowed to have endpoints with the same sign if the initial point has
078 * opposite sign function-wise).</p>
079 *
080 * @param f function to solve.
081 * @param min the lower bound for the interval.
082 * @param max the upper bound for the interval.
083 * @param initial the start value to use (must be set to min if no
084 * initial point is known).
085 * @return the value where the function is zero
086 * @throws MaxIterationsExceededException the maximum iteration count
087 * is exceeded
088 * @throws FunctionEvaluationException if an error occurs evaluating
089 * the function
090 * @throws IllegalArgumentException if initial is not between min and max
091 * (even if it <em>is</em> a root)
092 */
093 public double solve(final UnivariateRealFunction f,
094 final double min, final double max, final double initial)
095 throws MaxIterationsExceededException, FunctionEvaluationException {
096
097 clearResult();
098 verifySequence(min, initial, max);
099
100 // return the initial guess if it is good enough
101 double yInitial = f.value(initial);
102 if (Math.abs(yInitial) <= functionValueAccuracy) {
103 setResult(initial, 0);
104 return result;
105 }
106
107 // return the first endpoint if it is good enough
108 double yMin = f.value(min);
109 if (Math.abs(yMin) <= functionValueAccuracy) {
110 setResult(yMin, 0);
111 return result;
112 }
113
114 // reduce interval if min and initial bracket the root
115 if (yInitial * yMin < 0) {
116 return solve(f, min, yMin, initial, yInitial, min, yMin);
117 }
118
119 // return the second endpoint if it is good enough
120 double yMax = f.value(max);
121 if (Math.abs(yMax) <= functionValueAccuracy) {
122 setResult(yMax, 0);
123 return result;
124 }
125
126 // reduce interval if initial and max bracket the root
127 if (yInitial * yMax < 0) {
128 return solve(f, initial, yInitial, max, yMax, initial, yInitial);
129 }
130
131 // full Brent algorithm starting with provided initial guess
132 return solve(f, min, yMin, max, yMax, initial, yInitial);
133
134 }
135
136 /**
137 * Find a zero in the given interval.
138 * <p>
139 * Requires that the values of the function at the endpoints have opposite
140 * signs. An <code>IllegalArgumentException</code> is thrown if this is not
141 * the case.</p>
142 *
143 * @param f the function to solve
144 * @param min the lower bound for the interval.
145 * @param max the upper bound for the interval.
146 * @return the value where the function is zero
147 * @throws MaxIterationsExceededException if the maximum iteration count is exceeded
148 * @throws FunctionEvaluationException if an error occurs evaluating the
149 * function
150 * @throws IllegalArgumentException if min is not less than max or the
151 * signs of the values of the function at the endpoints are not opposites
152 */
153 public double solve(final UnivariateRealFunction f,
154 final double min, final double max)
155 throws MaxIterationsExceededException,
156 FunctionEvaluationException {
157
158 clearResult();
159 verifyInterval(min, max);
160
161 double ret = Double.NaN;
162
163 double yMin = f.value(min);
164 double yMax = f.value(max);
165
166 // Verify bracketing
167 double sign = yMin * yMax;
168 if (sign > 0) {
169 // check if either value is close to a zero
170 if (Math.abs(yMin) <= functionValueAccuracy) {
171 setResult(min, 0);
172 ret = min;
173 } else if (Math.abs(yMax) <= functionValueAccuracy) {
174 setResult(max, 0);
175 ret = max;
176 } else {
177 // neither value is close to zero and min and max do not bracket root.
178 throw MathRuntimeException.createIllegalArgumentException(
179 "function values at endpoints do not have different signs. " +
180 "Endpoints: [{0}, {1}], Values: [{2}, {3}]",
181 min, max, yMin, yMax);
182 }
183 } else if (sign < 0){
184 // solve using only the first endpoint as initial guess
185 ret = solve(f, min, yMin, max, yMax, min, yMin);
186 } else {
187 // either min or max is a root
188 if (yMin == 0.0) {
189 ret = min;
190 } else {
191 ret = max;
192 }
193 }
194
195 return ret;
196 }
197
198 /**
199 * Find a zero starting search according to the three provided points.
200 * @param f the function to solve
201 * @param x0 old approximation for the root
202 * @param y0 function value at the approximation for the root
203 * @param x1 last calculated approximation for the root
204 * @param y1 function value at the last calculated approximation
205 * for the root
206 * @param x2 bracket point (must be set to x0 if no bracket point is
207 * known, this will force starting with linear interpolation)
208 * @param y2 function value at the bracket point.
209 * @return the value where the function is zero
210 * @throws MaxIterationsExceededException if the maximum iteration count
211 * is exceeded
212 * @throws FunctionEvaluationException if an error occurs evaluating
213 * the function
214 */
215 private double solve(final UnivariateRealFunction f,
216 double x0, double y0,
217 double x1, double y1,
218 double x2, double y2)
219 throws MaxIterationsExceededException, FunctionEvaluationException {
220
221 double delta = x1 - x0;
222 double oldDelta = delta;
223
224 int i = 0;
225 while (i < maximalIterationCount) {
226 if (Math.abs(y2) < Math.abs(y1)) {
227 // use the bracket point if is better than last approximation
228 x0 = x1;
229 x1 = x2;
230 x2 = x0;
231 y0 = y1;
232 y1 = y2;
233 y2 = y0;
234 }
235 if (Math.abs(y1) <= functionValueAccuracy) {
236 // Avoid division by very small values. Assume
237 // the iteration has converged (the problem may
238 // still be ill conditioned)
239 setResult(x1, i);
240 return result;
241 }
242 double dx = (x2 - x1);
243 double tolerance =
244 Math.max(relativeAccuracy * Math.abs(x1), absoluteAccuracy);
245 if (Math.abs(dx) <= tolerance) {
246 setResult(x1, i);
247 return result;
248 }
249 if ((Math.abs(oldDelta) < tolerance) ||
250 (Math.abs(y0) <= Math.abs(y1))) {
251 // Force bisection.
252 delta = 0.5 * dx;
253 oldDelta = delta;
254 } else {
255 double r3 = y1 / y0;
256 double p;
257 double p1;
258 // the equality test (x0 == x2) is intentional,
259 // it is part of the original Brent's method,
260 // it should NOT be replaced by proximity test
261 if (x0 == x2) {
262 // Linear interpolation.
263 p = dx * r3;
264 p1 = 1.0 - r3;
265 } else {
266 // Inverse quadratic interpolation.
267 double r1 = y0 / y2;
268 double r2 = y1 / y2;
269 p = r3 * (dx * r1 * (r1 - r2) - (x1 - x0) * (r2 - 1.0));
270 p1 = (r1 - 1.0) * (r2 - 1.0) * (r3 - 1.0);
271 }
272 if (p > 0.0) {
273 p1 = -p1;
274 } else {
275 p = -p;
276 }
277 if (2.0 * p >= 1.5 * dx * p1 - Math.abs(tolerance * p1) ||
278 p >= Math.abs(0.5 * oldDelta * p1)) {
279 // Inverse quadratic interpolation gives a value
280 // in the wrong direction, or progress is slow.
281 // Fall back to bisection.
282 delta = 0.5 * dx;
283 oldDelta = delta;
284 } else {
285 oldDelta = delta;
286 delta = p / p1;
287 }
288 }
289 // Save old X1, Y1
290 x0 = x1;
291 y0 = y1;
292 // Compute new X1, Y1
293 if (Math.abs(delta) > tolerance) {
294 x1 = x1 + delta;
295 } else if (dx > 0.0) {
296 x1 = x1 + 0.5 * tolerance;
297 } else if (dx <= 0.0) {
298 x1 = x1 - 0.5 * tolerance;
299 }
300 y1 = f.value(x1);
301 if ((y1 > 0) == (y2 > 0)) {
302 x2 = x0;
303 y2 = y0;
304 delta = x1 - x0;
305 oldDelta = delta;
306 }
307 i++;
308 }
309 throw new MaxIterationsExceededException(maximalIterationCount);
310 }
311 }