    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    }