    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.MathRuntimeException;
    022    import org.apache.commons.math.MaxIterationsExceededException;
    023    import org.apache.commons.math.analysis.UnivariateRealFunction;
    024    import org.apache.commons.math.exception.util.LocalizedFormats;
    025    import org.apache.commons.math.util.FastMath;
    026
    027
    028    /**
    029     * Implements a modified version of the
    030     * <a href="http://mathworld.wolfram.com/SecantMethod.html">secant method</a>
    031     * for approximating a zero of a real univariate function.
    032     * <p>
    033     * The algorithm is modified to maintain bracketing of a root by successive
    034     * approximations. Because of forced bracketing, convergence may be slower than
    035     * the unrestricted secant algorithm. However, this implementation should in
    036     * general outperform the
    037     * <a href="http://mathworld.wolfram.com/MethodofFalsePosition.html">
    038     * regula falsi method.</a></p>
    039     * <p>
    040     * The function is assumed to be continuous but not necessarily smooth.</p>
    041     *
    042     * @version $Revision: 1070725 $ $Date: 2011-02-15 02:31:12 +0100 (mar. 15 f??vr. 2011) $
    043     */
    044    public class SecantSolver extends UnivariateRealSolverImpl {
    045
    046        /**
    047         * Construct a solver for the given function.
    048         * @param f function to solve.
    049         * @deprecated as of 2.0 the function to solve is passed as an argument
    050         * to the {@link #solve(UnivariateRealFunction, double, double)} or
    051         * {@link UnivariateRealSolverImpl#solve(UnivariateRealFunction, double, double, double)}
    052         * method.
    053         */
    054        @Deprecated
    055        public SecantSolver(UnivariateRealFunction f) {
    056            super(f, 100, 1E-6);
    057        }
    058
    059        /**
    060         * Construct a solver.
    061         * @deprecated in 2.2 (to be removed in 3.0).
    062         */
    063        @Deprecated
    064        public SecantSolver() {
    065            super(100, 1E-6);
    066        }
    067
    068        /** {@inheritDoc} */
    069        @Deprecated
    070        public double solve(final double min, final double max)
    071            throws ConvergenceException, FunctionEvaluationException {
    072            return solve(f, min, max);
    073        }
    074
    075        /** {@inheritDoc} */
    076        @Deprecated
    077        public double solve(final double min, final double max, final double initial)
    078            throws ConvergenceException, FunctionEvaluationException {
    079            return solve(f, min, max, initial);
    080        }
    081
    082        /**
    083         * Find a zero in the given interval.
    084         *
    085         * @param f the function to solve
    086         * @param min the lower bound for the interval
    087         * @param max the upper bound for the interval
    088         * @param initial the start value to use (ignored)
    089         * @param maxEval Maximum number of evaluations.
    090         * @return the value where the function is zero
    091         * @throws MaxIterationsExceededException if the maximum iteration count is exceeded
    092         * @throws FunctionEvaluationException if an error occurs evaluating the function
    093         * @throws IllegalArgumentException if min is not less than max or the
    094         * signs of the values of the function at the endpoints are not opposites
    095         */
    096        @Override
    097        public double solve(int maxEval, final UnivariateRealFunction f,
    098                            final double min, final double max, final double initial)
    099            throws MaxIterationsExceededException, FunctionEvaluationException {
    100            setMaximalIterationCount(maxEval);
    101            return solve(f, min, max, initial);
    102        }
    103
    104        /**
    105         * Find a zero in the given interval.
    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 (ignored)
    111         * @return the value where the function is zero
    112         * @throws MaxIterationsExceededException if the maximum iteration count is exceeded
    113         * @throws FunctionEvaluationException if an error occurs evaluating the function
    114         * @throws IllegalArgumentException if min is not less than max or the
    115         * signs of the values of the function at the endpoints are not opposites
    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            return solve(f, min, max);
    123        }
    124
    125        /**
    126         * Find a zero in the given interval.
    127         * @param f the function to solve
    128         * @param min the lower bound for the interval.
    129         * @param max the upper bound for the interval.
    130         * @param maxEval Maximum number of evaluations.
    131         * @return the value where the function is zero
    132         * @throws MaxIterationsExceededException  if the maximum iteration count is exceeded
    133         * @throws FunctionEvaluationException if an error occurs evaluating the function
    134         * @throws IllegalArgumentException if min is not less than max or the
    135         * signs of the values of the function at the endpoints are not opposites
    136         */
    137        @Override
    138        public double solve(int maxEval, final UnivariateRealFunction f,
    139                            final double min, final double max)
    140            throws MaxIterationsExceededException, FunctionEvaluationException {
    141            setMaximalIterationCount(maxEval);
    142            return solve(f, min, max);
    143        }
    144
    145        /**
    146         * Find a zero in the given interval.
    147         * @param f the function to solve
    148         * @param min the lower bound for the interval.
    149         * @param max the upper bound for the interval.
    150         * @return the value where the function is zero
    151         * @throws MaxIterationsExceededException  if the maximum iteration count is exceeded
    152         * @throws FunctionEvaluationException if an error occurs evaluating the function
    153         * @throws IllegalArgumentException if min is not less than max or the
    154         * signs of the values of the function at the endpoints are not opposites
    155         * @deprecated in 2.2 (to be removed in 3.0).
    156         */
    157        @Deprecated
    158        public double solve(final UnivariateRealFunction f,
    159                            final double min, final double max)
    160            throws MaxIterationsExceededException, FunctionEvaluationException {
    161
    162            clearResult();
    163            verifyInterval(min, max);
    164
    165            // Index 0 is the old approximation for the root.
    166            // Index 1 is the last calculated approximation  for the root.
    167            // Index 2 is a bracket for the root with respect to x0.
    168            // OldDelta is the length of the bracketing interval of the last
    169            // iteration.
    170            double x0 = min;
    171            double x1 = max;
    172            double y0 = f.value(x0);
    173            double y1 = f.value(x1);
    174
    175            // Verify bracketing
    176            if (y0 * y1 >= 0) {
    177                throw MathRuntimeException.createIllegalArgumentException(
    178                      LocalizedFormats.SAME_SIGN_AT_ENDPOINTS, min, max, y0, y1);
    179            }
    180
    181            double x2 = x0;
    182            double y2 = y0;
    183            double oldDelta = x2 - x1;
    184            int i = 0;
    185            while (i < maximalIterationCount) {
    186                if (FastMath.abs(y2) < FastMath.abs(y1)) {
    187                    x0 = x1;
    188                    x1 = x2;
    189                    x2 = x0;
    190                    y0 = y1;
    191                    y1 = y2;
    192                    y2 = y0;
    193                }
    194                if (FastMath.abs(y1) <= functionValueAccuracy) {
    195                    setResult(x1, i);
    196                    return result;
    197                }
    198                if (FastMath.abs(oldDelta) <
    199                    FastMath.max(relativeAccuracy * FastMath.abs(x1), absoluteAccuracy)) {
    200                    setResult(x1, i);
    201                    return result;
    202                }
    203                double delta;
    204                if (FastMath.abs(y1) > FastMath.abs(y0)) {
    205                    // Function value increased in last iteration. Force bisection.
    206                    delta = 0.5 * oldDelta;
    207                } else {
    208                    delta = (x0 - x1) / (1 - y0 / y1);
    209                    if (delta / oldDelta > 1) {
    210                        // New approximation falls outside bracket.
    211                        // Fall back to bisection.
    212                        delta = 0.5 * oldDelta;
    213                    }
    214                }
    215                x0 = x1;
    216                y0 = y1;
    217                x1 = x1 + delta;
    218                y1 = f.value(x1);
    219                if ((y1 > 0) == (y2 > 0)) {
    220                    // New bracket is (x0,x1).
    221                    x2 = x0;
    222                    y2 = y0;
    223                }
    224                oldDelta = x2 - x1;
    225                i++;
    226            }
    227            throw new MaxIterationsExceededException(maximalIterationCount);
    228        }
    229
    230    }