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.math3.analysis.solvers;
18  
19  import org.apache.commons.math3.util.FastMath;
20  import org.apache.commons.math3.exception.NumberIsTooLargeException;
21  import org.apache.commons.math3.exception.NoBracketingException;
22  import org.apache.commons.math3.exception.TooManyEvaluationsException;
23  
24  /**
25   * This class implements the <a href="http://mathworld.wolfram.com/MullersMethod.html">
26   * Muller's Method</a> for root finding of real univariate functions. For
27   * reference, see <b>Elementary Numerical Analysis</b>, ISBN 0070124477,
28   * chapter 3.
29   * <p>
30   * Muller's method applies to both real and complex functions, but here we
31   * restrict ourselves to real functions.
32   * This class differs from {@link MullerSolver} in the way it avoids complex
33   * operations.</p>
34   * Muller's original method would have function evaluation at complex point.
35   * Since our f(x) is real, we have to find ways to avoid that. Bracketing
36   * condition is one way to go: by requiring bracketing in every iteration,
37   * the newly computed approximation is guaranteed to be real.</p>
38   * <p>
39   * Normally Muller's method converges quadratically in the vicinity of a
40   * zero, however it may be very slow in regions far away from zeros. For
41   * example, f(x) = exp(x) - 1, min = -50, max = 100. In such case we use
42   * bisection as a safety backup if it performs very poorly.</p>
43   * <p>
44   * The formulas here use divided differences directly.</p>
45   *
46   * @since 1.2
47   * @see MullerSolver2
48   */
49  public class MullerSolver extends AbstractUnivariateSolver {
50  
51      /** Default absolute accuracy. */
52      private static final double DEFAULT_ABSOLUTE_ACCURACY = 1e-6;
53  
54      /**
55       * Construct a solver with default accuracy (1e-6).
56       */
57      public MullerSolver() {
58          this(DEFAULT_ABSOLUTE_ACCURACY);
59      }
60      /**
61       * Construct a solver.
62       *
63       * @param absoluteAccuracy Absolute accuracy.
64       */
65      public MullerSolver(double absoluteAccuracy) {
66          super(absoluteAccuracy);
67      }
68      /**
69       * Construct a solver.
70       *
71       * @param relativeAccuracy Relative accuracy.
72       * @param absoluteAccuracy Absolute accuracy.
73       */
74      public MullerSolver(double relativeAccuracy,
75                          double absoluteAccuracy) {
76          super(relativeAccuracy, absoluteAccuracy);
77      }
78  
79      /**
80       * {@inheritDoc}
81       */
82      @Override
83      protected double doSolve()
84          throws TooManyEvaluationsException,
85                 NumberIsTooLargeException,
86                 NoBracketingException {
87          final double min = getMin();
88          final double max = getMax();
89          final double initial = getStartValue();
90  
91          final double functionValueAccuracy = getFunctionValueAccuracy();
92  
93          verifySequence(min, initial, max);
94  
95          // check for zeros before verifying bracketing
96          final double fMin = computeObjectiveValue(min);
97          if (FastMath.abs(fMin) < functionValueAccuracy) {
98              return min;
99          }
100         final double fMax = computeObjectiveValue(max);
101         if (FastMath.abs(fMax) < functionValueAccuracy) {
102             return max;
103         }
104         final double fInitial = computeObjectiveValue(initial);
105         if (FastMath.abs(fInitial) <  functionValueAccuracy) {
106             return initial;
107         }
108 
109         verifyBracketing(min, max);
110 
111         if (isBracketing(min, initial)) {
112             return solve(min, initial, fMin, fInitial);
113         } else {
114             return solve(initial, max, fInitial, fMax);
115         }
116     }
117 
118     /**
119      * Find a real root in the given interval.
120      *
121      * @param min Lower bound for the interval.
122      * @param max Upper bound for the interval.
123      * @param fMin function value at the lower bound.
124      * @param fMax function value at the upper bound.
125      * @return the point at which the function value is zero.
126      * @throws TooManyEvaluationsException if the allowed number of calls to
127      * the function to be solved has been exhausted.
128      */
129     private double solve(double min, double max,
130                          double fMin, double fMax)
131         throws TooManyEvaluationsException {
132         final double relativeAccuracy = getRelativeAccuracy();
133         final double absoluteAccuracy = getAbsoluteAccuracy();
134         final double functionValueAccuracy = getFunctionValueAccuracy();
135 
136         // [x0, x2] is the bracketing interval in each iteration
137         // x1 is the last approximation and an interpolation point in (x0, x2)
138         // x is the new root approximation and new x1 for next round
139         // d01, d12, d012 are divided differences
140 
141         double x0 = min;
142         double y0 = fMin;
143         double x2 = max;
144         double y2 = fMax;
145         double x1 = 0.5 * (x0 + x2);
146         double y1 = computeObjectiveValue(x1);
147 
148         double oldx = Double.POSITIVE_INFINITY;
149         while (true) {
150             // Muller's method employs quadratic interpolation through
151             // x0, x1, x2 and x is the zero of the interpolating parabola.
152             // Due to bracketing condition, this parabola must have two
153             // real roots and we choose one in [x0, x2] to be x.
154             final double d01 = (y1 - y0) / (x1 - x0);
155             final double d12 = (y2 - y1) / (x2 - x1);
156             final double d012 = (d12 - d01) / (x2 - x0);
157             final double c1 = d01 + (x1 - x0) * d012;
158             final double delta = c1 * c1 - 4 * y1 * d012;
159             final double xplus = x1 + (-2.0 * y1) / (c1 + FastMath.sqrt(delta));
160             final double xminus = x1 + (-2.0 * y1) / (c1 - FastMath.sqrt(delta));
161             // xplus and xminus are two roots of parabola and at least
162             // one of them should lie in (x0, x2)
163             final double x = isSequence(x0, xplus, x2) ? xplus : xminus;
164             final double y = computeObjectiveValue(x);
165 
166             // check for convergence
167             final double tolerance = FastMath.max(relativeAccuracy * FastMath.abs(x), absoluteAccuracy);
168             if (FastMath.abs(x - oldx) <= tolerance ||
169                 FastMath.abs(y) <= functionValueAccuracy) {
170                 return x;
171             }
172 
173             // Bisect if convergence is too slow. Bisection would waste
174             // our calculation of x, hopefully it won't happen often.
175             // the real number equality test x == x1 is intentional and
176             // completes the proximity tests above it
177             boolean bisect = (x < x1 && (x1 - x0) > 0.95 * (x2 - x0)) ||
178                              (x > x1 && (x2 - x1) > 0.95 * (x2 - x0)) ||
179                              (x == x1);
180             // prepare the new bracketing interval for next iteration
181             if (!bisect) {
182                 x0 = x < x1 ? x0 : x1;
183                 y0 = x < x1 ? y0 : y1;
184                 x2 = x > x1 ? x2 : x1;
185                 y2 = x > x1 ? y2 : y1;
186                 x1 = x; y1 = y;
187                 oldx = x;
188             } else {
189                 double xm = 0.5 * (x0 + x2);
190                 double ym = computeObjectiveValue(xm);
191                 if (FastMath.signum(y0) + FastMath.signum(ym) == 0.0) {
192                     x2 = xm; y2 = ym;
193                 } else {
194                     x0 = xm; y0 = ym;
195                 }
196                 x1 = 0.5 * (x0 + x2);
197                 y1 = computeObjectiveValue(x1);
198                 oldx = Double.POSITIVE_INFINITY;
199             }
200         }
201     }
202 }