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 import org.apache.commons.math.ConvergenceException;
020 import org.apache.commons.math.FunctionEvaluationException;
021 import org.apache.commons.math.MaxIterationsExceededException;
022 import org.apache.commons.math.analysis.UnivariateRealFunction;
023 import org.apache.commons.math.util.MathUtils;
024
025 /**
026 * Implements the <a href="http://mathworld.wolfram.com/MullersMethod.html">
027 * Muller's Method</a> for root finding of real univariate functions. For
028 * reference, see <b>Elementary Numerical Analysis</b>, ISBN 0070124477,
029 * chapter 3.
030 * <p>
031 * Muller's method applies to both real and complex functions, but here we
032 * restrict ourselves to real functions. Methods solve() and solve2() find
033 * real zeros, using different ways to bypass complex arithmetics.</p>
034 *
035 * @version $Revision: 799857 $ $Date: 2009-08-01 09:07:12 -0400 (Sat, 01 Aug 2009) $
036 * @since 1.2
037 */
038 public class MullerSolver extends UnivariateRealSolverImpl {
039
040 /**
041 * Construct a solver for the given function.
042 *
043 * @param f function to solve
044 * @deprecated as of 2.0 the function to solve is passed as an argument
045 * to the {@link #solve(UnivariateRealFunction, double, double)} or
046 * {@link UnivariateRealSolverImpl#solve(UnivariateRealFunction, double, double, double)}
047 * method.
048 */
049 @Deprecated
050 public MullerSolver(UnivariateRealFunction f) {
051 super(f, 100, 1E-6);
052 }
053
054 /**
055 * Construct a solver.
056 */
057 public MullerSolver() {
058 super(100, 1E-6);
059 }
060
061 /** {@inheritDoc} */
062 @Deprecated
063 public double solve(final double min, final double max)
064 throws ConvergenceException, FunctionEvaluationException {
065 return solve(f, min, max);
066 }
067
068 /** {@inheritDoc} */
069 @Deprecated
070 public double solve(final double min, final double max, final double initial)
071 throws ConvergenceException, FunctionEvaluationException {
072 return solve(f, min, max, initial);
073 }
074
075 /**
076 * Find a real root in the given interval with initial value.
077 * <p>
078 * Requires bracketing condition.</p>
079 *
080 * @param f the 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
084 * @return the point at which the function value is zero
085 * @throws MaxIterationsExceededException if the maximum iteration count is exceeded
086 * or the solver detects convergence problems otherwise
087 * @throws FunctionEvaluationException if an error occurs evaluating the
088 * function
089 * @throws IllegalArgumentException if any parameters are invalid
090 */
091 public double solve(final UnivariateRealFunction f,
092 final double min, final double max, final double initial)
093 throws MaxIterationsExceededException, FunctionEvaluationException {
094
095 // check for zeros before verifying bracketing
096 if (f.value(min) == 0.0) { return min; }
097 if (f.value(max) == 0.0) { return max; }
098 if (f.value(initial) == 0.0) { return initial; }
099
100 verifyBracketing(min, max, f);
101 verifySequence(min, initial, max);
102 if (isBracketing(min, initial, f)) {
103 return solve(f, min, initial);
104 } else {
105 return solve(f, initial, max);
106 }
107 }
108
109 /**
110 * Find a real root in the given interval.
111 * <p>
112 * Original Muller's method would have function evaluation at complex point.
113 * Since our f(x) is real, we have to find ways to avoid that. Bracketing
114 * condition is one way to go: by requiring bracketing in every iteration,
115 * the newly computed approximation is guaranteed to be real.</p>
116 * <p>
117 * Normally Muller's method converges quadratically in the vicinity of a
118 * zero, however it may be very slow in regions far away from zeros. For
119 * example, f(x) = exp(x) - 1, min = -50, max = 100. In such case we use
120 * bisection as a safety backup if it performs very poorly.</p>
121 * <p>
122 * The formulas here use divided differences directly.</p>
123 *
124 * @param f the function to solve
125 * @param min the lower bound for the interval
126 * @param max the upper bound for the interval
127 * @return the point at which the function value is zero
128 * @throws MaxIterationsExceededException if the maximum iteration count is exceeded
129 * or the solver detects convergence problems otherwise
130 * @throws FunctionEvaluationException if an error occurs evaluating the
131 * function
132 * @throws IllegalArgumentException if any parameters are invalid
133 */
134 public double solve(final UnivariateRealFunction f,
135 final double min, final double max)
136 throws MaxIterationsExceededException, FunctionEvaluationException {
137
138 // [x0, x2] is the bracketing interval in each iteration
139 // x1 is the last approximation and an interpolation point in (x0, x2)
140 // x is the new root approximation and new x1 for next round
141 // d01, d12, d012 are divided differences
142 double x0, x1, x2, x, oldx, y0, y1, y2, y;
143 double d01, d12, d012, c1, delta, xplus, xminus, tolerance;
144
145 x0 = min; y0 = f.value(x0);
146 x2 = max; y2 = f.value(x2);
147 x1 = 0.5 * (x0 + x2); y1 = f.value(x1);
148
149 // check for zeros before verifying bracketing
150 if (y0 == 0.0) { return min; }
151 if (y2 == 0.0) { return max; }
152 verifyBracketing(min, max, f);
153
154 int i = 1;
155 oldx = Double.POSITIVE_INFINITY;
156 while (i <= maximalIterationCount) {
157 // Muller's method employs quadratic interpolation through
158 // x0, x1, x2 and x is the zero of the interpolating parabola.
159 // Due to bracketing condition, this parabola must have two
160 // real roots and we choose one in [x0, x2] to be x.
161 d01 = (y1 - y0) / (x1 - x0);
162 d12 = (y2 - y1) / (x2 - x1);
163 d012 = (d12 - d01) / (x2 - x0);
164 c1 = d01 + (x1 - x0) * d012;
165 delta = c1 * c1 - 4 * y1 * d012;
166 xplus = x1 + (-2.0 * y1) / (c1 + Math.sqrt(delta));
167 xminus = x1 + (-2.0 * y1) / (c1 - Math.sqrt(delta));
168 // xplus and xminus are two roots of parabola and at least
169 // one of them should lie in (x0, x2)
170 x = isSequence(x0, xplus, x2) ? xplus : xminus;
171 y = f.value(x);
172
173 // check for convergence
174 tolerance = Math.max(relativeAccuracy * Math.abs(x), absoluteAccuracy);
175 if (Math.abs(x - oldx) <= tolerance) {
176 setResult(x, i);
177 return result;
178 }
179 if (Math.abs(y) <= functionValueAccuracy) {
180 setResult(x, i);
181 return result;
182 }
183
184 // Bisect if convergence is too slow. Bisection would waste
185 // our calculation of x, hopefully it won't happen often.
186 // the real number equality test x == x1 is intentional and
187 // completes the proximity tests above it
188 boolean bisect = (x < x1 && (x1 - x0) > 0.95 * (x2 - x0)) ||
189 (x > x1 && (x2 - x1) > 0.95 * (x2 - x0)) ||
190 (x == x1);
191 // prepare the new bracketing interval for next iteration
192 if (!bisect) {
193 x0 = x < x1 ? x0 : x1;
194 y0 = x < x1 ? y0 : y1;
195 x2 = x > x1 ? x2 : x1;
196 y2 = x > x1 ? y2 : y1;
197 x1 = x; y1 = y;
198 oldx = x;
199 } else {
200 double xm = 0.5 * (x0 + x2);
201 double ym = f.value(xm);
202 if (MathUtils.sign(y0) + MathUtils.sign(ym) == 0.0) {
203 x2 = xm; y2 = ym;
204 } else {
205 x0 = xm; y0 = ym;
206 }
207 x1 = 0.5 * (x0 + x2);
208 y1 = f.value(x1);
209 oldx = Double.POSITIVE_INFINITY;
210 }
211 i++;
212 }
213 throw new MaxIterationsExceededException(maximalIterationCount);
214 }
215
216 /**
217 * Find a real root in the given interval.
218 * <p>
219 * solve2() differs from solve() in the way it avoids complex operations.
220 * Except for the initial [min, max], solve2() does not require bracketing
221 * condition, e.g. f(x0), f(x1), f(x2) can have the same sign. If complex
222 * number arises in the computation, we simply use its modulus as real
223 * approximation.</p>
224 * <p>
225 * Because the interval may not be bracketing, bisection alternative is
226 * not applicable here. However in practice our treatment usually works
227 * well, especially near real zeros where the imaginary part of complex
228 * approximation is often negligible.</p>
229 * <p>
230 * The formulas here do not use divided differences directly.</p>
231 *
232 * @param min the lower bound for the interval
233 * @param max the upper bound for the interval
234 * @return the point at which the function value is zero
235 * @throws MaxIterationsExceededException if the maximum iteration count is exceeded
236 * or the solver detects convergence problems otherwise
237 * @throws FunctionEvaluationException if an error occurs evaluating the
238 * function
239 * @throws IllegalArgumentException if any parameters are invalid
240 * @deprecated replaced by {@link #solve2(UnivariateRealFunction, double, double)}
241 * since 2.0
242 */
243 @Deprecated
244 public double solve2(final double min, final double max)
245 throws MaxIterationsExceededException, FunctionEvaluationException {
246 return solve2(f, min, max);
247 }
248
249 /**
250 * Find a real root in the given interval.
251 * <p>
252 * solve2() differs from solve() in the way it avoids complex operations.
253 * Except for the initial [min, max], solve2() does not require bracketing
254 * condition, e.g. f(x0), f(x1), f(x2) can have the same sign. If complex
255 * number arises in the computation, we simply use its modulus as real
256 * approximation.</p>
257 * <p>
258 * Because the interval may not be bracketing, bisection alternative is
259 * not applicable here. However in practice our treatment usually works
260 * well, especially near real zeros where the imaginary part of complex
261 * approximation is often negligible.</p>
262 * <p>
263 * The formulas here do not use divided differences directly.</p>
264 *
265 * @param f the function to solve
266 * @param min the lower bound for the interval
267 * @param max the upper bound for the interval
268 * @return the point at which the function value is zero
269 * @throws MaxIterationsExceededException if the maximum iteration count is exceeded
270 * or the solver detects convergence problems otherwise
271 * @throws FunctionEvaluationException if an error occurs evaluating the
272 * function
273 * @throws IllegalArgumentException if any parameters are invalid
274 */
275 public double solve2(final UnivariateRealFunction f,
276 final double min, final double max)
277 throws MaxIterationsExceededException, FunctionEvaluationException {
278
279 // x2 is the last root approximation
280 // x is the new approximation and new x2 for next round
281 // x0 < x1 < x2 does not hold here
282 double x0, x1, x2, x, oldx, y0, y1, y2, y;
283 double q, A, B, C, delta, denominator, tolerance;
284
285 x0 = min; y0 = f.value(x0);
286 x1 = max; y1 = f.value(x1);
287 x2 = 0.5 * (x0 + x1); y2 = f.value(x2);
288
289 // check for zeros before verifying bracketing
290 if (y0 == 0.0) { return min; }
291 if (y1 == 0.0) { return max; }
292 verifyBracketing(min, max, f);
293
294 int i = 1;
295 oldx = Double.POSITIVE_INFINITY;
296 while (i <= maximalIterationCount) {
297 // quadratic interpolation through x0, x1, x2
298 q = (x2 - x1) / (x1 - x0);
299 A = q * (y2 - (1 + q) * y1 + q * y0);
300 B = (2*q + 1) * y2 - (1 + q) * (1 + q) * y1 + q * q * y0;
301 C = (1 + q) * y2;
302 delta = B * B - 4 * A * C;
303 if (delta >= 0.0) {
304 // choose a denominator larger in magnitude
305 double dplus = B + Math.sqrt(delta);
306 double dminus = B - Math.sqrt(delta);
307 denominator = Math.abs(dplus) > Math.abs(dminus) ? dplus : dminus;
308 } else {
309 // take the modulus of (B +/- Math.sqrt(delta))
310 denominator = Math.sqrt(B * B - delta);
311 }
312 if (denominator != 0) {
313 x = x2 - 2.0 * C * (x2 - x1) / denominator;
314 // perturb x if it exactly coincides with x1 or x2
315 // the equality tests here are intentional
316 while (x == x1 || x == x2) {
317 x += absoluteAccuracy;
318 }
319 } else {
320 // extremely rare case, get a random number to skip it
321 x = min + Math.random() * (max - min);
322 oldx = Double.POSITIVE_INFINITY;
323 }
324 y = f.value(x);
325
326 // check for convergence
327 tolerance = Math.max(relativeAccuracy * Math.abs(x), absoluteAccuracy);
328 if (Math.abs(x - oldx) <= tolerance) {
329 setResult(x, i);
330 return result;
331 }
332 if (Math.abs(y) <= functionValueAccuracy) {
333 setResult(x, i);
334 return result;
335 }
336
337 // prepare the next iteration
338 x0 = x1; y0 = y1;
339 x1 = x2; y1 = y2;
340 x2 = x; y2 = y;
341 oldx = x;
342 i++;
343 }
344 throw new MaxIterationsExceededException(maximalIterationCount);
345 }
346 }