    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.exception.NoBracketingException;
    020    import org.apache.commons.math3.exception.NumberIsTooLargeException;
    021    import org.apache.commons.math3.exception.TooManyEvaluationsException;
    022    import org.apache.commons.math3.util.FastMath;
    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     * Except for the initial [min, max], it does not require bracketing
    035     * condition, e.g. f(x0), f(x1), f(x2) can have the same sign. If complex
    036     * number arises in the computation, we simply use its modulus as real
    037     * approximation.</p>
    038     * <p>
    039     * Because the interval may not be bracketing, bisection alternative is
    040     * not applicable here. However in practice our treatment usually works
    041     * well, especially near real zeroes where the imaginary part of complex
    042     * approximation is often negligible.</p>
    043     * <p>
    044     * The formulas here do not use divided differences directly.</p>
    045     *
    046     * @version $Id: MullerSolver2.java 1379560 2012-08-31 19:40:30Z erans $
    047     * @since 1.2
    048     * @see MullerSolver
    049     */
    050    public class MullerSolver2 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 MullerSolver2() {
    059            this(DEFAULT_ABSOLUTE_ACCURACY);
    060        }
    061        /**
    062         * Construct a solver.
    063         *
    064         * @param absoluteAccuracy Absolute accuracy.
    065         */
    066        public MullerSolver2(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 MullerSolver2(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
    091            verifyInterval(min, max);
    092
    093            final double relativeAccuracy = getRelativeAccuracy();
    094            final double absoluteAccuracy = getAbsoluteAccuracy();
    095            final double functionValueAccuracy = getFunctionValueAccuracy();
    096
    097            // x2 is the last root approximation
    098            // x is the new approximation and new x2 for next round
    099            // x0 < x1 < x2 does not hold here
    100
    101            double x0 = min;
    102            double y0 = computeObjectiveValue(x0);
    103            if (FastMath.abs(y0) < functionValueAccuracy) {
    104                return x0;
    105            }
    106            double x1 = max;
    107            double y1 = computeObjectiveValue(x1);
    108            if (FastMath.abs(y1) < functionValueAccuracy) {
    109                return x1;
    110            }
    111
    112            if(y0 * y1 > 0) {
    113                throw new NoBracketingException(x0, x1, y0, y1);
    114            }
    115
    116            double x2 = 0.5 * (x0 + x1);
    117            double y2 = computeObjectiveValue(x2);
    118
    119            double oldx = Double.POSITIVE_INFINITY;
    120            while (true) {
    121                // quadratic interpolation through x0, x1, x2
    122                final double q = (x2 - x1) / (x1 - x0);
    123                final double a = q * (y2 - (1 + q) * y1 + q * y0);
    124                final double b = (2 * q + 1) * y2 - (1 + q) * (1 + q) * y1 + q * q * y0;
    125                final double c = (1 + q) * y2;
    126                final double delta = b * b - 4 * a * c;
    127                double x;
    128                final double denominator;
    129                if (delta >= 0.0) {
    130                    // choose a denominator larger in magnitude
    131                    double dplus = b + FastMath.sqrt(delta);
    132                    double dminus = b - FastMath.sqrt(delta);
    133                    denominator = FastMath.abs(dplus) > FastMath.abs(dminus) ? dplus : dminus;
    134                } else {
    135                    // take the modulus of (B +/- FastMath.sqrt(delta))
    136                    denominator = FastMath.sqrt(b * b - delta);
    137                }
    138                if (denominator != 0) {
    139                    x = x2 - 2.0 * c * (x2 - x1) / denominator;
    140                    // perturb x if it exactly coincides with x1 or x2
    141                    // the equality tests here are intentional
    142                    while (x == x1 || x == x2) {
    143                        x += absoluteAccuracy;
    144                    }
    145                } else {
    146                    // extremely rare case, get a random number to skip it
    147                    x = min + FastMath.random() * (max - min);
    148                    oldx = Double.POSITIVE_INFINITY;
    149                }
    150                final double y = computeObjectiveValue(x);
    151
    152                // check for convergence
    153                final double tolerance = FastMath.max(relativeAccuracy * FastMath.abs(x), absoluteAccuracy);
    154                if (FastMath.abs(x - oldx) <= tolerance ||
    155                    FastMath.abs(y) <= functionValueAccuracy) {
    156                    return x;
    157                }
    158
    159                // prepare the next iteration
    160                x0 = x1;
    161                y0 = y1;
    162                x1 = x2;
    163                y1 = y2;
    164                x2 = x;
    165                y2 = y;
    166                oldx = x;
    167            }
    168        }
    169    }