1 /*
2 * Licensed to the Apache Software Foundation (ASF) under one or more
3 * contributor license agreements. See the NOTICE file distributed with
4 * this work for additional information regarding copyright ownership.
5 * The ASF licenses this file to You under the Apache License, Version 2.0
6 * (the "License"); you may not use this file except in compliance with
7 * the License. You may obtain a copy of the License at
8 *
9 * http://www.apache.org/licenses/LICENSE-2.0
10 *
11 * Unless required by applicable law or agreed to in writing, software
12 * distributed under the License is distributed on an "AS IS" BASIS,
13 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14 * See the License for the specific language governing permissions and
15 * limitations under the License.
16 */
17 package org.apache.commons.math.analysis.solvers;
18
19 import org.apache.commons.math.ConvergenceException;
20 import org.apache.commons.math.FunctionEvaluationException;
21 import org.apache.commons.math.MaxIterationsExceededException;
22 import org.apache.commons.math.analysis.UnivariateRealFunction;
23 import org.apache.commons.math.util.MathUtils;
24
25 /**
26 * Implements the <a href="http://mathworld.wolfram.com/MullersMethod.html">
27 * Muller's Method</a> for root finding of real univariate functions. For
28 * reference, see <b>Elementary Numerical Analysis</b>, ISBN 0070124477,
29 * chapter 3.
30 * <p>
31 * Muller's method applies to both real and complex functions, but here we
32 * restrict ourselves to real functions. Methods solve() and solve2() find
33 * real zeros, using different ways to bypass complex arithmetics.</p>
34 *
35 * @version $Revision: 825919 $ $Date: 2009-10-16 10:51:55 -0400 (Fri, 16 Oct 2009) $
36 * @since 1.2
37 */
38 public class MullerSolver extends UnivariateRealSolverImpl {
39
40 /**
41 * Construct a solver for the given function.
42 *
43 * @param f function to solve
44 * @deprecated as of 2.0 the function to solve is passed as an argument
45 * to the {@link #solve(UnivariateRealFunction, double, double)} or
46 * {@link UnivariateRealSolverImpl#solve(UnivariateRealFunction, double, double, double)}
47 * method.
48 */
49 @Deprecated
50 public MullerSolver(UnivariateRealFunction f) {
51 super(f, 100, 1E-6);
52 }
53
54 /**
55 * Construct a solver.
56 */
57 public MullerSolver() {
58 super(100, 1E-6);
59 }
60
61 /** {@inheritDoc} */
62 @Deprecated
63 public double solve(final double min, final double max)
64 throws ConvergenceException, FunctionEvaluationException {
65 return solve(f, min, max);
66 }
67
68 /** {@inheritDoc} */
69 @Deprecated
70 public double solve(final double min, final double max, final double initial)
71 throws ConvergenceException, FunctionEvaluationException {
72 return solve(f, min, max, initial);
73 }
74
75 /**
76 * Find a real root in the given interval with initial value.
77 * <p>
78 * Requires bracketing condition.</p>
79 *
80 * @param f the function to solve
81 * @param min the lower bound for the interval
82 * @param max the upper bound for the interval
83 * @param initial the start value to use
84 * @return the point at which the function value is zero
85 * @throws MaxIterationsExceededException if the maximum iteration count is exceeded
86 * or the solver detects convergence problems otherwise
87 * @throws FunctionEvaluationException if an error occurs evaluating the
88 * function
89 * @throws IllegalArgumentException if any parameters are invalid
90 */
91 public double solve(final UnivariateRealFunction f,
92 final double min, final double max, final double initial)
93 throws MaxIterationsExceededException, FunctionEvaluationException {
94
95 // check for zeros before verifying bracketing
96 if (f.value(min) == 0.0) { return min; }
97 if (f.value(max) == 0.0) { return max; }
98 if (f.value(initial) == 0.0) { return initial; }
99
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
143 double x0 = min;
144 double y0 = f.value(x0);
145 double x2 = max;
146 double y2 = f.value(x2);
147 double x1 = 0.5 * (x0 + x2);
148 double y1 = f.value(x1);
149
150 // check for zeros before verifying bracketing
151 if (y0 == 0.0) {
152 return min;
153 }
154 if (y2 == 0.0) {
155 return max;
156 }
157 verifyBracketing(min, max, f);
158
159 double oldx = Double.POSITIVE_INFINITY;
160 for (int i = 1; i <= maximalIterationCount; ++i) {
161 // Muller's method employs quadratic interpolation through
162 // x0, x1, x2 and x is the zero of the interpolating parabola.
163 // Due to bracketing condition, this parabola must have two
164 // real roots and we choose one in [x0, x2] to be x.
165 final double d01 = (y1 - y0) / (x1 - x0);
166 final double d12 = (y2 - y1) / (x2 - x1);
167 final double d012 = (d12 - d01) / (x2 - x0);
168 final double c1 = d01 + (x1 - x0) * d012;
169 final double delta = c1 * c1 - 4 * y1 * d012;
170 final double xplus = x1 + (-2.0 * y1) / (c1 + Math.sqrt(delta));
171 final double xminus = x1 + (-2.0 * y1) / (c1 - Math.sqrt(delta));
172 // xplus and xminus are two roots of parabola and at least
173 // one of them should lie in (x0, x2)
174 final double x = isSequence(x0, xplus, x2) ? xplus : xminus;
175 final double y = f.value(x);
176
177 // check for convergence
178 final double tolerance = Math.max(relativeAccuracy * Math.abs(x), absoluteAccuracy);
179 if (Math.abs(x - oldx) <= tolerance) {
180 setResult(x, i);
181 return result;
182 }
183 if (Math.abs(y) <= functionValueAccuracy) {
184 setResult(x, i);
185 return result;
186 }
187
188 // Bisect if convergence is too slow. Bisection would waste
189 // our calculation of x, hopefully it won't happen often.
190 // the real number equality test x == x1 is intentional and
191 // completes the proximity tests above it
192 boolean bisect = (x < x1 && (x1 - x0) > 0.95 * (x2 - x0)) ||
193 (x > x1 && (x2 - x1) > 0.95 * (x2 - x0)) ||
194 (x == x1);
195 // prepare the new bracketing interval for next iteration
196 if (!bisect) {
197 x0 = x < x1 ? x0 : x1;
198 y0 = x < x1 ? y0 : y1;
199 x2 = x > x1 ? x2 : x1;
200 y2 = x > x1 ? y2 : y1;
201 x1 = x; y1 = y;
202 oldx = x;
203 } else {
204 double xm = 0.5 * (x0 + x2);
205 double ym = f.value(xm);
206 if (MathUtils.sign(y0) + MathUtils.sign(ym) == 0.0) {
207 x2 = xm; y2 = ym;
208 } else {
209 x0 = xm; y0 = ym;
210 }
211 x1 = 0.5 * (x0 + x2);
212 y1 = f.value(x1);
213 oldx = Double.POSITIVE_INFINITY;
214 }
215 }
216 throw new MaxIterationsExceededException(maximalIterationCount);
217 }
218
219 /**
220 * Find a real root in the given interval.
221 * <p>
222 * solve2() differs from solve() in the way it avoids complex operations.
223 * Except for the initial [min, max], solve2() does not require bracketing
224 * condition, e.g. f(x0), f(x1), f(x2) can have the same sign. If complex
225 * number arises in the computation, we simply use its modulus as real
226 * approximation.</p>
227 * <p>
228 * Because the interval may not be bracketing, bisection alternative is
229 * not applicable here. However in practice our treatment usually works
230 * well, especially near real zeros where the imaginary part of complex
231 * approximation is often negligible.</p>
232 * <p>
233 * The formulas here do not use divided differences directly.</p>
234 *
235 * @param min the lower bound for the interval
236 * @param max the upper bound for the interval
237 * @return the point at which the function value is zero
238 * @throws MaxIterationsExceededException if the maximum iteration count is exceeded
239 * or the solver detects convergence problems otherwise
240 * @throws FunctionEvaluationException if an error occurs evaluating the
241 * function
242 * @throws IllegalArgumentException if any parameters are invalid
243 * @deprecated replaced by {@link #solve2(UnivariateRealFunction, double, double)}
244 * since 2.0
245 */
246 @Deprecated
247 public double solve2(final double min, final double max)
248 throws MaxIterationsExceededException, FunctionEvaluationException {
249 return solve2(f, min, max);
250 }
251
252 /**
253 * Find a real root in the given interval.
254 * <p>
255 * solve2() differs from solve() in the way it avoids complex operations.
256 * Except for the initial [min, max], solve2() does not require bracketing
257 * condition, e.g. f(x0), f(x1), f(x2) can have the same sign. If complex
258 * number arises in the computation, we simply use its modulus as real
259 * approximation.</p>
260 * <p>
261 * Because the interval may not be bracketing, bisection alternative is
262 * not applicable here. However in practice our treatment usually works
263 * well, especially near real zeros where the imaginary part of complex
264 * approximation is often negligible.</p>
265 * <p>
266 * The formulas here do not use divided differences directly.</p>
267 *
268 * @param f the function to solve
269 * @param min the lower bound for the interval
270 * @param max the upper bound for the interval
271 * @return the point at which the function value is zero
272 * @throws MaxIterationsExceededException if the maximum iteration count is exceeded
273 * or the solver detects convergence problems otherwise
274 * @throws FunctionEvaluationException if an error occurs evaluating the
275 * function
276 * @throws IllegalArgumentException if any parameters are invalid
277 */
278 public double solve2(final UnivariateRealFunction f,
279 final double min, final double max)
280 throws MaxIterationsExceededException, FunctionEvaluationException {
281
282 // x2 is the last root approximation
283 // x is the new approximation and new x2 for next round
284 // x0 < x1 < x2 does not hold here
285
286 double x0 = min;
287 double y0 = f.value(x0);
288 double x1 = max;
289 double y1 = f.value(x1);
290 double x2 = 0.5 * (x0 + x1);
291 double y2 = f.value(x2);
292
293 // check for zeros before verifying bracketing
294 if (y0 == 0.0) { return min; }
295 if (y1 == 0.0) { return max; }
296 verifyBracketing(min, max, f);
297
298 double oldx = Double.POSITIVE_INFINITY;
299 for (int i = 1; i <= maximalIterationCount; ++i) {
300 // quadratic interpolation through x0, x1, x2
301 final double q = (x2 - x1) / (x1 - x0);
302 final double a = q * (y2 - (1 + q) * y1 + q * y0);
303 final double b = (2 * q + 1) * y2 - (1 + q) * (1 + q) * y1 + q * q * y0;
304 final double c = (1 + q) * y2;
305 final double delta = b * b - 4 * a * c;
306 double x;
307 final double denominator;
308 if (delta >= 0.0) {
309 // choose a denominator larger in magnitude
310 double dplus = b + Math.sqrt(delta);
311 double dminus = b - Math.sqrt(delta);
312 denominator = Math.abs(dplus) > Math.abs(dminus) ? dplus : dminus;
313 } else {
314 // take the modulus of (B +/- Math.sqrt(delta))
315 denominator = Math.sqrt(b * b - delta);
316 }
317 if (denominator != 0) {
318 x = x2 - 2.0 * c * (x2 - x1) / denominator;
319 // perturb x if it exactly coincides with x1 or x2
320 // the equality tests here are intentional
321 while (x == x1 || x == x2) {
322 x += absoluteAccuracy;
323 }
324 } else {
325 // extremely rare case, get a random number to skip it
326 x = min + Math.random() * (max - min);
327 oldx = Double.POSITIVE_INFINITY;
328 }
329 final double y = f.value(x);
330
331 // check for convergence
332 final double tolerance = Math.max(relativeAccuracy * Math.abs(x), absoluteAccuracy);
333 if (Math.abs(x - oldx) <= tolerance) {
334 setResult(x, i);
335 return result;
336 }
337 if (Math.abs(y) <= functionValueAccuracy) {
338 setResult(x, i);
339 return result;
340 }
341
342 // prepare the next iteration
343 x0 = x1;
344 y0 = y1;
345 x1 = x2;
346 y1 = y2;
347 x2 = x;
348 y2 = y;
349 oldx = x;
350 }
351 throw new MaxIterationsExceededException(maximalIterationCount);
352 }
353 }