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    }