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    }