    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.math.analysis.solvers;
    018
    019    import org.apache.commons.math.ConvergenceException;
    020    import org.apache.commons.math.FunctionEvaluationException;
    021    import org.apache.commons.math.MaxIterationsExceededException;
    022    import org.apache.commons.math.analysis.UnivariateRealFunction;
    023    import org.apache.commons.math.util.FastMath;
    024    import org.apache.commons.math.util.MathUtils;
    025
    026    /**
    027     * Implements the <a href="http://mathworld.wolfram.com/MullersMethod.html">
    028     * Muller's Method</a> for root finding of real univariate functions. For
    029     * reference, see <b>Elementary Numerical Analysis</b>, ISBN 0070124477,
    030     * chapter 3.
    031     * <p>
    032     * Muller's method applies to both real and complex functions, but here we
    033     * restrict ourselves to real functions. Methods solve() and solve2() find
    034     * real zeros, using different ways to bypass complex arithmetics.</p>
    035     *
    036     * @version $Revision: 1070725 $ $Date: 2011-02-15 02:31:12 +0100 (mar. 15 f??vr. 2011) $
    037     * @since 1.2
    038     */
    039    public class MullerSolver extends UnivariateRealSolverImpl {
    040
    041        /**
    042         * Construct a solver for the given function.
    043         *
    044         * @param f function to solve
    045         * @deprecated as of 2.0 the function to solve is passed as an argument
    046         * to the {@link #solve(UnivariateRealFunction, double, double)} or
    047         * {@link UnivariateRealSolverImpl#solve(UnivariateRealFunction, double, double, double)}
    048         * method.
    049         */
    050        @Deprecated
    051        public MullerSolver(UnivariateRealFunction f) {
    052            super(f, 100, 1E-6);
    053        }
    054
    055        /**
    056         * Construct a solver.
    057         * @deprecated in 2.2 (to be removed in 3.0).
    058         */
    059        @Deprecated
    060        public MullerSolver() {
    061            super(100, 1E-6);
    062        }
    063
    064        /** {@inheritDoc} */
    065        @Deprecated
    066        public double solve(final double min, final double max)
    067            throws ConvergenceException, FunctionEvaluationException {
    068            return solve(f, min, max);
    069        }
    070
    071        /** {@inheritDoc} */
    072        @Deprecated
    073        public double solve(final double min, final double max, final double initial)
    074            throws ConvergenceException, FunctionEvaluationException {
    075            return solve(f, min, max, initial);
    076        }
    077
    078        /**
    079         * Find a real root in the given interval with initial value.
    080         * <p>
    081         * Requires bracketing condition.</p>
    082         *
    083         * @param f the function to solve
    084         * @param min the lower bound for the interval
    085         * @param max the upper bound for the interval
    086         * @param initial the start value to use
    087         * @param maxEval Maximum number of evaluations.
    088         * @return the point at which the function value is zero
    089         * @throws MaxIterationsExceededException if the maximum iteration count is exceeded
    090         * or the solver detects convergence problems otherwise
    091         * @throws FunctionEvaluationException if an error occurs evaluating the function
    092         * @throws IllegalArgumentException if any parameters are invalid
    093         */
    094        @Override
    095        public double solve(int maxEval, final UnivariateRealFunction f,
    096                            final double min, final double max, final double initial)
    097            throws MaxIterationsExceededException, FunctionEvaluationException {
    098            setMaximalIterationCount(maxEval);
    099            return solve(f, min, max, initial);
    100        }
    101
    102        /**
    103         * Find a real root in the given interval with initial value.
    104         * <p>
    105         * Requires bracketing condition.</p>
    106         *
    107         * @param f the function to solve
    108         * @param min the lower bound for the interval
    109         * @param max the upper bound for the interval
    110         * @param initial the start value to use
    111         * @return the point at which the function value is zero
    112         * @throws MaxIterationsExceededException if the maximum iteration count is exceeded
    113         * or the solver detects convergence problems otherwise
    114         * @throws FunctionEvaluationException if an error occurs evaluating the function
    115         * @throws IllegalArgumentException if any parameters are invalid
    116         * @deprecated in 2.2 (to be removed in 3.0).
    117         */
    118        @Deprecated
    119        public double solve(final UnivariateRealFunction f,
    120                            final double min, final double max, final double initial)
    121            throws MaxIterationsExceededException, FunctionEvaluationException {
    122
    123            // check for zeros before verifying bracketing
    124            if (f.value(min) == 0.0) { return min; }
    125            if (f.value(max) == 0.0) { return max; }
    126            if (f.value(initial) == 0.0) { return initial; }
    127
    128            verifyBracketing(min, max, f);
    129            verifySequence(min, initial, max);
    130            if (isBracketing(min, initial, f)) {
    131                return solve(f, min, initial);
    132            } else {
    133                return solve(f, initial, max);
    134            }
    135        }
    136
    137        /**
    138         * Find a real root in the given interval.
    139         * <p>
    140         * Original Muller's method would have function evaluation at complex point.
    141         * Since our f(x) is real, we have to find ways to avoid that. Bracketing
    142         * condition is one way to go: by requiring bracketing in every iteration,
    143         * the newly computed approximation is guaranteed to be real.</p>
    144         * <p>
    145         * Normally Muller's method converges quadratically in the vicinity of a
    146         * zero, however it may be very slow in regions far away from zeros. For
    147         * example, f(x) = exp(x) - 1, min = -50, max = 100. In such case we use
    148         * bisection as a safety backup if it performs very poorly.</p>
    149         * <p>
    150         * The formulas here use divided differences directly.</p>
    151         *
    152         * @param f the function to solve
    153         * @param min the lower bound for the interval
    154         * @param max the upper bound for the interval
    155         * @param maxEval Maximum number of evaluations.
    156         * @return the point at which the function value is zero
    157         * @throws MaxIterationsExceededException if the maximum iteration count is exceeded
    158         * or the solver detects convergence problems otherwise
    159         * @throws FunctionEvaluationException if an error occurs evaluating the function
    160         * @throws IllegalArgumentException if any parameters are invalid
    161         */
    162        @Override
    163        public double solve(int maxEval, final UnivariateRealFunction f,
    164                            final double min, final double max)
    165            throws MaxIterationsExceededException, FunctionEvaluationException {
    166            setMaximalIterationCount(maxEval);
    167            return solve(f, min, max);
    168        }
    169
    170        /**
    171         * Find a real root in the given interval.
    172         * <p>
    173         * Original Muller's method would have function evaluation at complex point.
    174         * Since our f(x) is real, we have to find ways to avoid that. Bracketing
    175         * condition is one way to go: by requiring bracketing in every iteration,
    176         * the newly computed approximation is guaranteed to be real.</p>
    177         * <p>
    178         * Normally Muller's method converges quadratically in the vicinity of a
    179         * zero, however it may be very slow in regions far away from zeros. For
    180         * example, f(x) = exp(x) - 1, min = -50, max = 100. In such case we use
    181         * bisection as a safety backup if it performs very poorly.</p>
    182         * <p>
    183         * The formulas here use divided differences directly.</p>
    184         *
    185         * @param f the function to solve
    186         * @param min the lower bound for the interval
    187         * @param max the upper bound for the interval
    188         * @return the point at which the function value is zero
    189         * @throws MaxIterationsExceededException if the maximum iteration count is exceeded
    190         * or the solver detects convergence problems otherwise
    191         * @throws FunctionEvaluationException if an error occurs evaluating the function
    192         * @throws IllegalArgumentException if any parameters are invalid
    193         * @deprecated in 2.2 (to be removed in 3.0).
    194         */
    195        @Deprecated
    196        public double solve(final UnivariateRealFunction f,
    197                            final double min, final double max)
    198            throws MaxIterationsExceededException, FunctionEvaluationException {
    199
    200            // [x0, x2] is the bracketing interval in each iteration
    201            // x1 is the last approximation and an interpolation point in (x0, x2)
    202            // x is the new root approximation and new x1 for next round
    203            // d01, d12, d012 are divided differences
    204
    205            double x0 = min;
    206            double y0 = f.value(x0);
    207            double x2 = max;
    208            double y2 = f.value(x2);
    209            double x1 = 0.5 * (x0 + x2);
    210            double y1 = f.value(x1);
    211
    212            // check for zeros before verifying bracketing
    213            if (y0 == 0.0) {
    214                return min;
    215            }
    216            if (y2 == 0.0) {
    217                return max;
    218            }
    219            verifyBracketing(min, max, f);
    220
    221            double oldx = Double.POSITIVE_INFINITY;
    222            for (int i = 1; i <= maximalIterationCount; ++i) {
    223                // Muller's method employs quadratic interpolation through
    224                // x0, x1, x2 and x is the zero of the interpolating parabola.
    225                // Due to bracketing condition, this parabola must have two
    226                // real roots and we choose one in [x0, x2] to be x.
    227                final double d01 = (y1 - y0) / (x1 - x0);
    228                final double d12 = (y2 - y1) / (x2 - x1);
    229                final double d012 = (d12 - d01) / (x2 - x0);
    230                final double c1 = d01 + (x1 - x0) * d012;
    231                final double delta = c1 * c1 - 4 * y1 * d012;
    232                final double xplus = x1 + (-2.0 * y1) / (c1 + FastMath.sqrt(delta));
    233                final double xminus = x1 + (-2.0 * y1) / (c1 - FastMath.sqrt(delta));
    234                // xplus and xminus are two roots of parabola and at least
    235                // one of them should lie in (x0, x2)
    236                final double x = isSequence(x0, xplus, x2) ? xplus : xminus;
    237                final double y = f.value(x);
    238
    239                // check for convergence
    240                final double tolerance = FastMath.max(relativeAccuracy * FastMath.abs(x), absoluteAccuracy);
    241                if (FastMath.abs(x - oldx) <= tolerance) {
    242                    setResult(x, i);
    243                    return result;
    244                }
    245                if (FastMath.abs(y) <= functionValueAccuracy) {
    246                    setResult(x, i);
    247                    return result;
    248                }
    249
    250                // Bisect if convergence is too slow. Bisection would waste
    251                // our calculation of x, hopefully it won't happen often.
    252                // the real number equality test x == x1 is intentional and
    253                // completes the proximity tests above it
    254                boolean bisect = (x < x1 && (x1 - x0) > 0.95 * (x2 - x0)) ||
    255                                 (x > x1 && (x2 - x1) > 0.95 * (x2 - x0)) ||
    256                                 (x == x1);
    257                // prepare the new bracketing interval for next iteration
    258                if (!bisect) {
    259                    x0 = x < x1 ? x0 : x1;
    260                    y0 = x < x1 ? y0 : y1;
    261                    x2 = x > x1 ? x2 : x1;
    262                    y2 = x > x1 ? y2 : y1;
    263                    x1 = x; y1 = y;
    264                    oldx = x;
    265                } else {
    266                    double xm = 0.5 * (x0 + x2);
    267                    double ym = f.value(xm);
    268                    if (MathUtils.sign(y0) + MathUtils.sign(ym) == 0.0) {
    269                        x2 = xm; y2 = ym;
    270                    } else {
    271                        x0 = xm; y0 = ym;
    272                    }
    273                    x1 = 0.5 * (x0 + x2);
    274                    y1 = f.value(x1);
    275                    oldx = Double.POSITIVE_INFINITY;
    276                }
    277            }
    278            throw new MaxIterationsExceededException(maximalIterationCount);
    279        }
    280
    281        /**
    282         * Find a real root in the given interval.
    283         * <p>
    284         * solve2() differs from solve() in the way it avoids complex operations.
    285         * Except for the initial [min, max], solve2() does not require bracketing
    286         * condition, e.g. f(x0), f(x1), f(x2) can have the same sign. If complex
    287         * number arises in the computation, we simply use its modulus as real
    288         * approximation.</p>
    289         * <p>
    290         * Because the interval may not be bracketing, bisection alternative is
    291         * not applicable here. However in practice our treatment usually works
    292         * well, especially near real zeros where the imaginary part of complex
    293         * approximation is often negligible.</p>
    294         * <p>
    295         * The formulas here do not use divided differences directly.</p>
    296         *
    297         * @param min the lower bound for the interval
    298         * @param max the upper bound for the interval
    299         * @return the point at which the function value is zero
    300         * @throws MaxIterationsExceededException if the maximum iteration count is exceeded
    301         * or the solver detects convergence problems otherwise
    302         * @throws FunctionEvaluationException if an error occurs evaluating the function
    303         * @throws IllegalArgumentException if any parameters are invalid
    304         * @deprecated replaced by {@link #solve2(UnivariateRealFunction, double, double)}
    305         * since 2.0
    306         */
    307        @Deprecated
    308        public double solve2(final double min, final double max)
    309            throws MaxIterationsExceededException, FunctionEvaluationException {
    310            return solve2(f, min, max);
    311        }
    312
    313        /**
    314         * Find a real root in the given interval.
    315         * <p>
    316         * solve2() differs from solve() in the way it avoids complex operations.
    317         * Except for the initial [min, max], solve2() does not require bracketing
    318         * condition, e.g. f(x0), f(x1), f(x2) can have the same sign. If complex
    319         * number arises in the computation, we simply use its modulus as real
    320         * approximation.</p>
    321         * <p>
    322         * Because the interval may not be bracketing, bisection alternative is
    323         * not applicable here. However in practice our treatment usually works
    324         * well, especially near real zeros where the imaginary part of complex
    325         * approximation is often negligible.</p>
    326         * <p>
    327         * The formulas here do not use divided differences directly.</p>
    328         *
    329         * @param f the function to solve
    330         * @param min the lower bound for the interval
    331         * @param max the upper bound for the interval
    332         * @return the point at which the function value is zero
    333         * @throws MaxIterationsExceededException if the maximum iteration count is exceeded
    334         * or the solver detects convergence problems otherwise
    335         * @throws FunctionEvaluationException if an error occurs evaluating the function
    336         * @throws IllegalArgumentException if any parameters are invalid
    337         * @deprecated in 2.2 (to be removed in 3.0).
    338         */
    339        @Deprecated
    340        public double solve2(final UnivariateRealFunction f,
    341                             final double min, final double max)
    342            throws MaxIterationsExceededException, FunctionEvaluationException {
    343
    344            // x2 is the last root approximation
    345            // x is the new approximation and new x2 for next round
    346            // x0 < x1 < x2 does not hold here
    347
    348            double x0 = min;
    349            double y0 = f.value(x0);
    350            double x1 = max;
    351            double y1 = f.value(x1);
    352            double x2 = 0.5 * (x0 + x1);
    353            double y2 = f.value(x2);
    354
    355            // check for zeros before verifying bracketing
    356            if (y0 == 0.0) { return min; }
    357            if (y1 == 0.0) { return max; }
    358            verifyBracketing(min, max, f);
    359
    360            double oldx = Double.POSITIVE_INFINITY;
    361            for (int i = 1; i <= maximalIterationCount; ++i) {
    362                // quadratic interpolation through x0, x1, x2
    363                final double q = (x2 - x1) / (x1 - x0);
    364                final double a = q * (y2 - (1 + q) * y1 + q * y0);
    365                final double b = (2 * q + 1) * y2 - (1 + q) * (1 + q) * y1 + q * q * y0;
    366                final double c = (1 + q) * y2;
    367                final double delta = b * b - 4 * a * c;
    368                double x;
    369                final double denominator;
    370                if (delta >= 0.0) {
    371                    // choose a denominator larger in magnitude
    372                    double dplus = b + FastMath.sqrt(delta);
    373                    double dminus = b - FastMath.sqrt(delta);
    374                    denominator = FastMath.abs(dplus) > FastMath.abs(dminus) ? dplus : dminus;
    375                } else {
    376                    // take the modulus of (B +/- FastMath.sqrt(delta))
    377                    denominator = FastMath.sqrt(b * b - delta);
    378                }
    379                if (denominator != 0) {
    380                    x = x2 - 2.0 * c * (x2 - x1) / denominator;
    381                    // perturb x if it exactly coincides with x1 or x2
    382                    // the equality tests here are intentional
    383                    while (x == x1 || x == x2) {
    384                        x += absoluteAccuracy;
    385                    }
    386                } else {
    387                    // extremely rare case, get a random number to skip it
    388                    x = min + FastMath.random() * (max - min);
    389                    oldx = Double.POSITIVE_INFINITY;
    390                }
    391                final double y = f.value(x);
    392
    393                // check for convergence
    394                final double tolerance = FastMath.max(relativeAccuracy * FastMath.abs(x), absoluteAccuracy);
    395                if (FastMath.abs(x - oldx) <= tolerance) {
    396                    setResult(x, i);
    397                    return result;
    398                }
    399                if (FastMath.abs(y) <= functionValueAccuracy) {
    400                    setResult(x, i);
    401                    return result;
    402                }
    403
    404                // prepare the next iteration
    405                x0 = x1;
    406                y0 = y1;
    407                x1 = x2;
    408                y1 = y2;
    409                x2 = x;
    410                y2 = y;
    411                oldx = x;
    412            }
    413            throw new MaxIterationsExceededException(maximalIterationCount);
    414        }
    415    }