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 */
017package org.apache.commons.math3.analysis.solvers;
018
019import org.apache.commons.math3.exception.NoBracketingException;
020import org.apache.commons.math3.exception.NumberIsTooLargeException;
021import org.apache.commons.math3.exception.TooManyEvaluationsException;
022import 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><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 * @since 1.2
047 * @see MullerSolver2
048 */
049public class MullerSolver extends AbstractUnivariateSolver {
050
051    /** Default absolute accuracy. */
052    private static final double DEFAULT_ABSOLUTE_ACCURACY = 1e-6;
053
054    /**
055     * Construct a solver with default accuracy (1e-6).
056     */
057    public MullerSolver() {
058        this(DEFAULT_ABSOLUTE_ACCURACY);
059    }
060    /**
061     * Construct a solver.
062     *
063     * @param absoluteAccuracy Absolute accuracy.
064     */
065    public MullerSolver(double absoluteAccuracy) {
066        super(absoluteAccuracy);
067    }
068    /**
069     * Construct a solver.
070     *
071     * @param relativeAccuracy Relative accuracy.
072     * @param absoluteAccuracy Absolute accuracy.
073     */
074    public MullerSolver(double relativeAccuracy,
075                        double absoluteAccuracy) {
076        super(relativeAccuracy, absoluteAccuracy);
077    }
078
079    /**
080     * {@inheritDoc}
081     */
082    @Override
083    protected double doSolve()
084        throws TooManyEvaluationsException,
085               NumberIsTooLargeException,
086               NoBracketingException {
087        final double min = getMin();
088        final double max = getMax();
089        final double initial = getStartValue();
090
091        final double functionValueAccuracy = getFunctionValueAccuracy();
092
093        verifySequence(min, initial, max);
094
095        // check for zeros before verifying bracketing
096        final double fMin = computeObjectiveValue(min);
097        if (FastMath.abs(fMin) < functionValueAccuracy) {
098            return min;
099        }
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}