    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
    018    package org.apache.commons.math3.analysis.solvers;
    019
    020    import org.apache.commons.math3.util.FastMath;
    021    import org.apache.commons.math3.exception.NoBracketingException;
    022    import org.apache.commons.math3.exception.TooManyEvaluationsException;
    023
    024    /**
    025     * Implements the <em>Secant</em> method for root-finding (approximating a
    026     * zero of a univariate real function). The solution that is maintained is
    027     * not bracketed, and as such convergence is not guaranteed.
    028     *
    029     * <p>Implementation based on the following article: M. Dowell and P. Jarratt,
    030     * <em>A modified regula falsi method for computing the root of an
    031     * equation</em>, BIT Numerical Mathematics, volume 11, number 2,
    032     * pages 168-174, Springer, 1971.</p>
    033     *
    034     * <p>Note that since release 3.0 this class implements the actual
    035     * <em>Secant</em> algorithm, and not a modified one. As such, the 3.0 version
    036     * is not backwards compatible with previous versions. To use an algorithm
    037     * similar to the pre-3.0 releases, use the
    038     * {@link IllinoisSolver <em>Illinois</em>} algorithm or the
    039     * {@link PegasusSolver <em>Pegasus</em>} algorithm.</p>
    040     *
    041     * @version $Id: SecantSolver.java 1379560 2012-08-31 19:40:30Z erans $
    042     */
    043    public class SecantSolver extends AbstractUnivariateSolver {
    044
    045        /** Default absolute accuracy. */
    046        protected static final double DEFAULT_ABSOLUTE_ACCURACY = 1e-6;
    047
    048        /** Construct a solver with default accuracy (1e-6). */
    049        public SecantSolver() {
    050            super(DEFAULT_ABSOLUTE_ACCURACY);
    051        }
    052
    053        /**
    054         * Construct a solver.
    055         *
    056         * @param absoluteAccuracy absolute accuracy
    057         */
    058        public SecantSolver(final double absoluteAccuracy) {
    059            super(absoluteAccuracy);
    060        }
    061
    062        /**
    063         * Construct a solver.
    064         *
    065         * @param relativeAccuracy relative accuracy
    066         * @param absoluteAccuracy absolute accuracy
    067         */
    068        public SecantSolver(final double relativeAccuracy,
    069                            final double absoluteAccuracy) {
    070            super(relativeAccuracy, absoluteAccuracy);
    071        }
    072
    073        /** {@inheritDoc} */
    074        @Override
    075        protected final double doSolve()
    076            throws TooManyEvaluationsException,
    077                   NoBracketingException {
    078            // Get initial solution
    079            double x0 = getMin();
    080            double x1 = getMax();
    081            double f0 = computeObjectiveValue(x0);
    082            double f1 = computeObjectiveValue(x1);
    083
    084            // If one of the bounds is the exact root, return it. Since these are
    085            // not under-approximations or over-approximations, we can return them
    086            // regardless of the allowed solutions.
    087            if (f0 == 0.0) {
    088                return x0;
    089            }
    090            if (f1 == 0.0) {
    091                return x1;
    092            }
    093
    094            // Verify bracketing of initial solution.
    095            verifyBracketing(x0, x1);
    096
    097            // Get accuracies.
    098            final double ftol = getFunctionValueAccuracy();
    099            final double atol = getAbsoluteAccuracy();
    100            final double rtol = getRelativeAccuracy();
    101
    102            // Keep finding better approximations.
    103            while (true) {
    104                // Calculate the next approximation.
    105                final double x = x1 - ((f1 * (x1 - x0)) / (f1 - f0));
    106                final double fx = computeObjectiveValue(x);
    107
    108                // If the new approximation is the exact root, return it. Since
    109                // this is not an under-approximation or an over-approximation,
    110                // we can return it regardless of the allowed solutions.
    111                if (fx == 0.0) {
    112                    return x;
    113                }
    114
    115                // Update the bounds with the new approximation.
    116                x0 = x1;
    117                f0 = f1;
    118                x1 = x;
    119                f1 = fx;
    120
    121                // If the function value of the last approximation is too small,
    122                // given the function value accuracy, then we can't get closer to
    123                // the root than we already are.
    124                if (FastMath.abs(f1) <= ftol) {
    125                    return x1;
    126                }
    127
    128                // If the current interval is within the given accuracies, we
    129                // are satisfied with the current approximation.
    130                if (FastMath.abs(x1 - x0) < FastMath.max(rtol * FastMath.abs(x1), atol)) {
    131                    return x1;
    132                }
    133            }
    134        }
    135
    136    }