View Javadoc

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: 799857 $ $Date: 2009-08-01 09:07:12 -0400 (Sat, 01 Aug 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         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 }