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