    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/RiddersMethod.html">
    028     * Ridders' Method</a> for root finding of real univariate functions. For
    029     * reference, see C. Ridders, <i>A new algorithm for computing a single root
    030     * of a real continuous function </i>, IEEE Transactions on Circuits and
    031     * Systems, 26 (1979), 979 - 980.
    032     * <p>
    033     * The function should be continuous but not necessarily smooth.</p>
    034     *
    035     * @version $Revision: 1070725 $ $Date: 2011-02-15 02:31:12 +0100 (mar. 15 f??vr. 2011) $
    036     * @since 1.2
    037     */
    038    public class RiddersSolver 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 RiddersSolver(UnivariateRealFunction f) {
    051            super(f, 100, 1E-6);
    052        }
    053
    054        /**
    055         * Construct a solver.
    056         * @deprecated in 2.2
    057         */
    058        @Deprecated
    059        public RiddersSolver() {
    060            super(100, 1E-6);
    061        }
    062
    063        /** {@inheritDoc} */
    064        @Deprecated
    065        public double solve(final double min, final double max)
    066            throws ConvergenceException, FunctionEvaluationException {
    067            return solve(f, min, max);
    068        }
    069
    070        /** {@inheritDoc} */
    071        @Deprecated
    072        public double solve(final double min, final double max, final double initial)
    073            throws ConvergenceException, FunctionEvaluationException {
    074            return solve(f, min, max, initial);
    075        }
    076
    077        /**
    078         * Find a root in the given interval with initial value.
    079         * <p>
    080         * Requires bracketing condition.</p>
    081         *
    082         * @param f the function to solve
    083         * @param min the lower bound for the interval
    084         * @param max the upper bound for the interval
    085         * @param initial the start value to use
    086         * @param maxEval Maximum number of evaluations.
    087         * @return the point at which the function value is zero
    088         * @throws MaxIterationsExceededException if the maximum iteration count is exceeded
    089         * @throws FunctionEvaluationException if an error occurs evaluating the function
    090         * @throws IllegalArgumentException if any parameters are invalid
    091         */
    092        @Override
    093        public double solve(int maxEval, final UnivariateRealFunction f,
    094                            final double min, final double max, final double initial)
    095            throws MaxIterationsExceededException, FunctionEvaluationException {
    096            setMaximalIterationCount(maxEval);
    097            return solve(f, min, max, initial);
    098        }
    099
    100        /**
    101         * Find a root in the given interval with initial value.
    102         * <p>
    103         * Requires bracketing condition.</p>
    104         *
    105         * @param f the function to solve
    106         * @param min the lower bound for the interval
    107         * @param max the upper bound for the interval
    108         * @param initial the start value to use
    109         * @return the point at which the function value is zero
    110         * @throws MaxIterationsExceededException if the maximum iteration count is exceeded
    111         * @throws FunctionEvaluationException if an error occurs evaluating the function
    112         * @throws IllegalArgumentException if any parameters are invalid
    113         * @deprecated in 2.2 (to be removed in 3.0).
    114         */
    115        @Deprecated
    116        public double solve(final UnivariateRealFunction f,
    117                            final double min, final double max, final double initial)
    118            throws MaxIterationsExceededException, FunctionEvaluationException {
    119
    120            // check for zeros before verifying bracketing
    121            if (f.value(min) == 0.0) { return min; }
    122            if (f.value(max) == 0.0) { return max; }
    123            if (f.value(initial) == 0.0) { return initial; }
    124
    125            verifyBracketing(min, max, f);
    126            verifySequence(min, initial, max);
    127            if (isBracketing(min, initial, f)) {
    128                return solve(f, min, initial);
    129            } else {
    130                return solve(f, initial, max);
    131            }
    132        }
    133
    134        /**
    135         * Find a root in the given interval.
    136         * <p>
    137         * Requires bracketing condition.</p>
    138         *
    139         * @param f the function to solve
    140         * @param min the lower bound for the interval
    141         * @param max the upper bound for the interval
    142         * @param maxEval Maximum number of evaluations.
    143         * @return the point at which the function value is zero
    144         * @throws MaxIterationsExceededException if the maximum iteration count is exceeded
    145         * @throws FunctionEvaluationException if an error occurs evaluating the function
    146         * @throws IllegalArgumentException if any parameters are invalid
    147         */
    148        @Override
    149        public double solve(int maxEval, final UnivariateRealFunction f,
    150                            final double min, final double max)
    151            throws MaxIterationsExceededException, FunctionEvaluationException {
    152            setMaximalIterationCount(maxEval);
    153            return solve(f, min, max);
    154        }
    155
    156        /**
    157         * Find a root in the given interval.
    158         * <p>
    159         * Requires bracketing condition.</p>
    160         *
    161         * @param f the function to solve
    162         * @param min the lower bound for the interval
    163         * @param max the upper bound for the interval
    164         * @return the point at which the function value is zero
    165         * @throws MaxIterationsExceededException if the maximum iteration count is exceeded
    166         * @throws FunctionEvaluationException if an error occurs evaluating the function
    167         * @throws IllegalArgumentException if any parameters are invalid
    168         * @deprecated in 2.2 (to be removed in 3.0).
    169         */
    170        @Deprecated
    171        public double solve(final UnivariateRealFunction f,
    172                            final double min, final double max)
    173            throws MaxIterationsExceededException, FunctionEvaluationException {
    174
    175            // [x1, x2] is the bracketing interval in each iteration
    176            // x3 is the midpoint of [x1, x2]
    177            // x is the new root approximation and an endpoint of the new interval
    178            double x1 = min;
    179            double y1 = f.value(x1);
    180            double x2 = max;
    181            double y2 = f.value(x2);
    182
    183            // check for zeros before verifying bracketing
    184            if (y1 == 0.0) {
    185                return min;
    186            }
    187            if (y2 == 0.0) {
    188                return max;
    189            }
    190            verifyBracketing(min, max, f);
    191
    192            int i = 1;
    193            double oldx = Double.POSITIVE_INFINITY;
    194            while (i <= maximalIterationCount) {
    195                // calculate the new root approximation
    196                final double x3 = 0.5 * (x1 + x2);
    197                final double y3 = f.value(x3);
    198                if (FastMath.abs(y3) <= functionValueAccuracy) {
    199                    setResult(x3, i);
    200                    return result;
    201                }
    202                final double delta = 1 - (y1 * y2) / (y3 * y3);  // delta > 1 due to bracketing
    203                final double correction = (MathUtils.sign(y2) * MathUtils.sign(y3)) *
    204                                          (x3 - x1) / FastMath.sqrt(delta);
    205                final double x = x3 - correction;                // correction != 0
    206                final double y = f.value(x);
    207
    208                // check for convergence
    209                final double tolerance = FastMath.max(relativeAccuracy * FastMath.abs(x), absoluteAccuracy);
    210                if (FastMath.abs(x - oldx) <= tolerance) {
    211                    setResult(x, i);
    212                    return result;
    213                }
    214                if (FastMath.abs(y) <= functionValueAccuracy) {
    215                    setResult(x, i);
    216                    return result;
    217                }
    218
    219                // prepare the new interval for next iteration
    220                // Ridders' method guarantees x1 < x < x2
    221                if (correction > 0.0) {             // x1 < x < x3
    222                    if (MathUtils.sign(y1) + MathUtils.sign(y) == 0.0) {
    223                        x2 = x;
    224                        y2 = y;
    225                    } else {
    226                        x1 = x;
    227                        x2 = x3;
    228                        y1 = y;
    229                        y2 = y3;
    230                    }
    231                } else {                            // x3 < x < x2
    232                    if (MathUtils.sign(y2) + MathUtils.sign(y) == 0.0) {
    233                        x1 = x;
    234                        y1 = y;
    235                    } else {
    236                        x1 = x3;
    237                        x2 = x;
    238                        y1 = y3;
    239                        y2 = y;
    240                    }
    241                }
    242                oldx = x;
    243                i++;
    244            }
    245            throw new MaxIterationsExceededException(maximalIterationCount);
    246        }
    247    }