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.math3.analysis.solvers;
018    
019    import org.apache.commons.math3.util.FastMath;
020    import org.apache.commons.math3.exception.NumberIsTooLargeException;
021    import org.apache.commons.math3.exception.NoBracketingException;
022    import org.apache.commons.math3.exception.TooManyEvaluationsException;
023    
024    /**
025     * This class implements the <a href="http://mathworld.wolfram.com/MullersMethod.html">
026     * Muller's Method</a> for root finding of real univariate functions. For
027     * reference, see <b>Elementary Numerical Analysis</b>, ISBN 0070124477,
028     * chapter 3.
029     * <p>
030     * Muller's method applies to both real and complex functions, but here we
031     * restrict ourselves to real functions.
032     * This class differs from {@link MullerSolver} in the way it avoids complex
033     * operations.</p>
034     * Muller's original method would have function evaluation at complex point.
035     * Since our f(x) is real, we have to find ways to avoid that. Bracketing
036     * condition is one way to go: by requiring bracketing in every iteration,
037     * the newly computed approximation is guaranteed to be real.</p>
038     * <p>
039     * Normally Muller's method converges quadratically in the vicinity of a
040     * zero, however it may be very slow in regions far away from zeros. For
041     * example, f(x) = exp(x) - 1, min = -50, max = 100. In such case we use
042     * bisection as a safety backup if it performs very poorly.</p>
043     * <p>
044     * The formulas here use divided differences directly.</p>
045     *
046     * @version $Id: MullerSolver.java 1391927 2012-09-30 00:03:30Z erans $
047     * @since 1.2
048     * @see MullerSolver2
049     */
050    public class MullerSolver extends AbstractUnivariateSolver {
051    
052        /** Default absolute accuracy. */
053        private static final double DEFAULT_ABSOLUTE_ACCURACY = 1e-6;
054    
055        /**
056         * Construct a solver with default accuracy (1e-6).
057         */
058        public MullerSolver() {
059            this(DEFAULT_ABSOLUTE_ACCURACY);
060        }
061        /**
062         * Construct a solver.
063         *
064         * @param absoluteAccuracy Absolute accuracy.
065         */
066        public MullerSolver(double absoluteAccuracy) {
067            super(absoluteAccuracy);
068        }
069        /**
070         * Construct a solver.
071         *
072         * @param relativeAccuracy Relative accuracy.
073         * @param absoluteAccuracy Absolute accuracy.
074         */
075        public MullerSolver(double relativeAccuracy,
076                            double absoluteAccuracy) {
077            super(relativeAccuracy, absoluteAccuracy);
078        }
079    
080        /**
081         * {@inheritDoc}
082         */
083        @Override
084        protected double doSolve()
085            throws TooManyEvaluationsException,
086                   NumberIsTooLargeException,
087                   NoBracketingException {
088            final double min = getMin();
089            final double max = getMax();
090            final double initial = getStartValue();
091    
092            final double functionValueAccuracy = getFunctionValueAccuracy();
093    
094            verifySequence(min, initial, max);
095    
096            // check for zeros before verifying bracketing
097            final double fMin = computeObjectiveValue(min);
098            if (FastMath.abs(fMin) < functionValueAccuracy) {
099                return min;
100            }
101            final double fMax = computeObjectiveValue(max);
102            if (FastMath.abs(fMax) < functionValueAccuracy) {
103                return max;
104            }
105            final double fInitial = computeObjectiveValue(initial);
106            if (FastMath.abs(fInitial) <  functionValueAccuracy) {
107                return initial;
108            }
109    
110            verifyBracketing(min, max);
111    
112            if (isBracketing(min, initial)) {
113                return solve(min, initial, fMin, fInitial);
114            } else {
115                return solve(initial, max, fInitial, fMax);
116            }
117        }
118    
119        /**
120         * Find a real root in the given interval.
121         *
122         * @param min Lower bound for the interval.
123         * @param max Upper bound for the interval.
124         * @param fMin function value at the lower bound.
125         * @param fMax function value at the upper bound.
126         * @return the point at which the function value is zero.
127         * @throws TooManyEvaluationsException if the allowed number of calls to
128         * the function to be solved has been exhausted.
129         */
130        private double solve(double min, double max,
131                             double fMin, double fMax)
132            throws TooManyEvaluationsException {
133            final double relativeAccuracy = getRelativeAccuracy();
134            final double absoluteAccuracy = getAbsoluteAccuracy();
135            final double functionValueAccuracy = getFunctionValueAccuracy();
136    
137            // [x0, x2] is the bracketing interval in each iteration
138            // x1 is the last approximation and an interpolation point in (x0, x2)
139            // x is the new root approximation and new x1 for next round
140            // d01, d12, d012 are divided differences
141    
142            double x0 = min;
143            double y0 = fMin;
144            double x2 = max;
145            double y2 = fMax;
146            double x1 = 0.5 * (x0 + x2);
147            double y1 = computeObjectiveValue(x1);
148    
149            double oldx = Double.POSITIVE_INFINITY;
150            while (true) {
151                // Muller's method employs quadratic interpolation through
152                // x0, x1, x2 and x is the zero of the interpolating parabola.
153                // Due to bracketing condition, this parabola must have two
154                // real roots and we choose one in [x0, x2] to be x.
155                final double d01 = (y1 - y0) / (x1 - x0);
156                final double d12 = (y2 - y1) / (x2 - x1);
157                final double d012 = (d12 - d01) / (x2 - x0);
158                final double c1 = d01 + (x1 - x0) * d012;
159                final double delta = c1 * c1 - 4 * y1 * d012;
160                final double xplus = x1 + (-2.0 * y1) / (c1 + FastMath.sqrt(delta));
161                final double xminus = x1 + (-2.0 * y1) / (c1 - FastMath.sqrt(delta));
162                // xplus and xminus are two roots of parabola and at least
163                // one of them should lie in (x0, x2)
164                final double x = isSequence(x0, xplus, x2) ? xplus : xminus;
165                final double y = computeObjectiveValue(x);
166    
167                // check for convergence
168                final double tolerance = FastMath.max(relativeAccuracy * FastMath.abs(x), absoluteAccuracy);
169                if (FastMath.abs(x - oldx) <= tolerance ||
170                    FastMath.abs(y) <= functionValueAccuracy) {
171                    return x;
172                }
173    
174                // Bisect if convergence is too slow. Bisection would waste
175                // our calculation of x, hopefully it won't happen often.
176                // the real number equality test x == x1 is intentional and
177                // completes the proximity tests above it
178                boolean bisect = (x < x1 && (x1 - x0) > 0.95 * (x2 - x0)) ||
179                                 (x > x1 && (x2 - x1) > 0.95 * (x2 - x0)) ||
180                                 (x == x1);
181                // prepare the new bracketing interval for next iteration
182                if (!bisect) {
183                    x0 = x < x1 ? x0 : x1;
184                    y0 = x < x1 ? y0 : y1;
185                    x2 = x > x1 ? x2 : x1;
186                    y2 = x > x1 ? y2 : y1;
187                    x1 = x; y1 = y;
188                    oldx = x;
189                } else {
190                    double xm = 0.5 * (x0 + x2);
191                    double ym = computeObjectiveValue(xm);
192                    if (FastMath.signum(y0) + FastMath.signum(ym) == 0.0) {
193                        x2 = xm; y2 = ym;
194                    } else {
195                        x0 = xm; y0 = ym;
196                    }
197                    x1 = 0.5 * (x0 + x2);
198                    y1 = computeObjectiveValue(x1);
199                    oldx = Double.POSITIVE_INFINITY;
200                }
201            }
202        }
203    }