001: /*
002: * Copyright 2003-2005 The Apache Software Foundation.
003: *
004: * Licensed under the Apache License, Version 2.0 (the "License");
005: * you may not use this file except in compliance with the License.
006: * You may obtain a copy of the License at
007: *
008: * http://www.apache.org/licenses/LICENSE-2.0
009: *
010: * Unless required by applicable law or agreed to in writing, software
011: * distributed under the License is distributed on an "AS IS" BASIS,
012: * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
013: * See the License for the specific language governing permissions and
014: * limitations under the License.
015: */
016: package org.apache.commons.math.analysis;
017:
018: import org.apache.commons.math.ConvergenceException;
019: import org.apache.commons.math.FunctionEvaluationException;
020:
021: /**
022: * Implements the <a href="http://mathworld.wolfram.com/BrentsMethod.html">
023: * Brent algorithm</a> for finding zeros of real univariate functions.
024: * <p>
025: * The function should be continuous but not necessarily smooth.
026: *
027: * @version $Revision: 348519 $ $Date: 2005-11-23 12:12:18 -0700 (Wed, 23 Nov 2005) $
028: */
029: public class BrentSolver extends UnivariateRealSolverImpl {
030:
031: /** Serializable version identifier */
032: private static final long serialVersionUID = 3350616277306882875L;
033:
034: /**
035: * Construct a solver for the given function.
036: *
037: * @param f function to solve.
038: */
039: public BrentSolver(UnivariateRealFunction f) {
040: super (f, 100, 1E-6);
041: }
042:
043: /**
044: * Find a zero in the given interval.
045: * <p>
046: * Throws <code>ConvergenceException</code> if the values of the function
047: * at the endpoints of the interval have the same sign.
048: *
049: * @param min the lower bound for the interval.
050: * @param max the upper bound for the interval.
051: * @param initial the start value to use (ignored).
052: * @return the value where the function is zero
053: * @throws ConvergenceException the maximum iteration count is exceeded
054: * @throws FunctionEvaluationException if an error occurs evaluating
055: * the function
056: * @throws IllegalArgumentException if initial is not between min and max
057: */
058: public double solve(double min, double max, double initial)
059: throws ConvergenceException, FunctionEvaluationException {
060:
061: return solve(min, max);
062: }
063:
064: /**
065: * Find a zero in the given interval.
066: * <p>
067: * Requires that the values of the function at the endpoints have opposite
068: * signs. An <code>IllegalArgumentException</code> is thrown if this is not
069: * the case.
070: *
071: * @param min the lower bound for the interval.
072: * @param max the upper bound for the interval.
073: * @return the value where the function is zero
074: * @throws ConvergenceException if the maximum iteration count is exceeded
075: * @throws FunctionEvaluationException if an error occurs evaluating the
076: * function
077: * @throws IllegalArgumentException if min is not less than max or the
078: * signs of the values of the function at the endpoints are not opposites
079: */
080: public double solve(double min, double max)
081: throws ConvergenceException, FunctionEvaluationException {
082:
083: clearResult();
084: verifyInterval(min, max);
085:
086: // Index 0 is the old approximation for the root.
087: // Index 1 is the last calculated approximation for the root.
088: // Index 2 is a bracket for the root with respect to x1.
089: double x0 = min;
090: double x1 = max;
091: double y0;
092: double y1;
093: y0 = f.value(x0);
094: y1 = f.value(x1);
095:
096: // Verify bracketing
097: if (y0 * y1 >= 0) {
098: throw new IllegalArgumentException(
099: "Function values at endpoints do not have different signs."
100: + " Endpoints: [" + min + "," + max + "]"
101: + " Values: [" + y0 + "," + y1 + "]");
102: }
103:
104: double x2 = x0;
105: double y2 = y0;
106: double delta = x1 - x0;
107: double oldDelta = delta;
108:
109: int i = 0;
110: while (i < maximalIterationCount) {
111: if (Math.abs(y2) < Math.abs(y1)) {
112: x0 = x1;
113: x1 = x2;
114: x2 = x0;
115: y0 = y1;
116: y1 = y2;
117: y2 = y0;
118: }
119: if (Math.abs(y1) <= functionValueAccuracy) {
120: // Avoid division by very small values. Assume
121: // the iteration has converged (the problem may
122: // still be ill conditioned)
123: setResult(x1, i);
124: return result;
125: }
126: double dx = (x2 - x1);
127: double tolerance = Math.max(
128: relativeAccuracy * Math.abs(x1), absoluteAccuracy);
129: if (Math.abs(dx) <= tolerance) {
130: setResult(x1, i);
131: return result;
132: }
133: if ((Math.abs(oldDelta) < tolerance)
134: || (Math.abs(y0) <= Math.abs(y1))) {
135: // Force bisection.
136: delta = 0.5 * dx;
137: oldDelta = delta;
138: } else {
139: double r3 = y1 / y0;
140: double p;
141: double p1;
142: if (x0 == x2) {
143: // Linear interpolation.
144: p = dx * r3;
145: p1 = 1.0 - r3;
146: } else {
147: // Inverse quadratic interpolation.
148: double r1 = y0 / y2;
149: double r2 = y1 / y2;
150: p = r3
151: * (dx * r1 * (r1 - r2) - (x1 - x0)
152: * (r2 - 1.0));
153: p1 = (r1 - 1.0) * (r2 - 1.0) * (r3 - 1.0);
154: }
155: if (p > 0.0) {
156: p1 = -p1;
157: } else {
158: p = -p;
159: }
160: if (2.0 * p >= 1.5 * dx * p1 - Math.abs(tolerance * p1)
161: || p >= Math.abs(0.5 * oldDelta * p1)) {
162: // Inverse quadratic interpolation gives a value
163: // in the wrong direction, or progress is slow.
164: // Fall back to bisection.
165: delta = 0.5 * dx;
166: oldDelta = delta;
167: } else {
168: oldDelta = delta;
169: delta = p / p1;
170: }
171: }
172: // Save old X1, Y1
173: x0 = x1;
174: y0 = y1;
175: // Compute new X1, Y1
176: if (Math.abs(delta) > tolerance) {
177: x1 = x1 + delta;
178: } else if (dx > 0.0) {
179: x1 = x1 + 0.5 * tolerance;
180: } else if (dx <= 0.0) {
181: x1 = x1 - 0.5 * tolerance;
182: }
183: y1 = f.value(x1);
184: if ((y1 > 0) == (y2 > 0)) {
185: x2 = x0;
186: y2 = y0;
187: delta = x1 - x0;
188: oldDelta = delta;
189: }
190: i++;
191: }
192: throw new ConvergenceException(
193: "Maximum number of iterations exceeded.");
194: }
195: }
|