    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.analysis.UnivariateFunction;
    022    import org.apache.commons.math3.exception.ConvergenceException;
    023    import org.apache.commons.math3.exception.MathInternalError;
    024
    025    /**
    026     * Base class for all bracketing <em>Secant</em>-based methods for root-finding
    027     * (approximating a zero of a univariate real function).
    028     *
    029     * <p>Implementation of the {@link RegulaFalsiSolver <em>Regula Falsi</em>} and
    030     * {@link IllinoisSolver <em>Illinois</em>} methods is based on the
    031     * following article: M. Dowell and P. Jarratt,
    032     * <em>A modified regula falsi method for computing the root of an
    033     * equation</em>, BIT Numerical Mathematics, volume 11, number 2,
    034     * pages 168-174, Springer, 1971.</p>
    035     *
    036     * <p>Implementation of the {@link PegasusSolver <em>Pegasus</em>} method is
    037     * based on the following article: M. Dowell and P. Jarratt,
    038     * <em>The "Pegasus" method for computing the root of an equation</em>,
    039     * BIT Numerical Mathematics, volume 12, number 4, pages 503-508, Springer,
    040     * 1972.</p>
    041     *
    042     * <p>The {@link SecantSolver <em>Secant</em>} method is <em>not</em> a
    043     * bracketing method, so it is not implemented here. It has a separate
    044     * implementation.</p>
    045     *
    046     * @since 3.0
    047     * @version $Id: BaseSecantSolver.java 1379560 2012-08-31 19:40:30Z erans $
    048     */
    049    public abstract class BaseSecantSolver
    050        extends AbstractUnivariateSolver
    051        implements BracketedUnivariateSolver<UnivariateFunction> {
    052
    053        /** Default absolute accuracy. */
    054        protected static final double DEFAULT_ABSOLUTE_ACCURACY = 1e-6;
    055
    056        /** The kinds of solutions that the algorithm may accept. */
    057        private AllowedSolution allowed;
    058
    059        /** The <em>Secant</em>-based root-finding method to use. */
    060        private final Method method;
    061
    062        /**
    063         * Construct a solver.
    064         *
    065         * @param absoluteAccuracy Absolute accuracy.
    066         * @param method <em>Secant</em>-based root-finding method to use.
    067         */
    068        protected BaseSecantSolver(final double absoluteAccuracy, final Method method) {
    069            super(absoluteAccuracy);
    070            this.allowed = AllowedSolution.ANY_SIDE;
    071            this.method = method;
    072        }
    073
    074        /**
    075         * Construct a solver.
    076         *
    077         * @param relativeAccuracy Relative accuracy.
    078         * @param absoluteAccuracy Absolute accuracy.
    079         * @param method <em>Secant</em>-based root-finding method to use.
    080         */
    081        protected BaseSecantSolver(final double relativeAccuracy,
    082                                   final double absoluteAccuracy,
    083                                   final Method method) {
    084            super(relativeAccuracy, absoluteAccuracy);
    085            this.allowed = AllowedSolution.ANY_SIDE;
    086            this.method = method;
    087        }
    088
    089        /**
    090         * Construct a solver.
    091         *
    092         * @param relativeAccuracy Maximum relative error.
    093         * @param absoluteAccuracy Maximum absolute error.
    094         * @param functionValueAccuracy Maximum function value error.
    095         * @param method <em>Secant</em>-based root-finding method to use
    096         */
    097        protected BaseSecantSolver(final double relativeAccuracy,
    098                                   final double absoluteAccuracy,
    099                                   final double functionValueAccuracy,
    100                                   final Method method) {
    101            super(relativeAccuracy, absoluteAccuracy, functionValueAccuracy);
    102            this.allowed = AllowedSolution.ANY_SIDE;
    103            this.method = method;
    104        }
    105
    106        /** {@inheritDoc} */
    107        public double solve(final int maxEval, final UnivariateFunction f,
    108                            final double min, final double max,
    109                            final AllowedSolution allowedSolution) {
    110            return solve(maxEval, f, min, max, min + 0.5 * (max - min), allowedSolution);
    111        }
    112
    113        /** {@inheritDoc} */
    114        public double solve(final int maxEval, final UnivariateFunction f,
    115                            final double min, final double max, final double startValue,
    116                            final AllowedSolution allowedSolution) {
    117            this.allowed = allowedSolution;
    118            return super.solve(maxEval, f, min, max, startValue);
    119        }
    120
    121        /** {@inheritDoc} */
    122        @Override
    123        public double solve(final int maxEval, final UnivariateFunction f,
    124                            final double min, final double max, final double startValue) {
    125            return solve(maxEval, f, min, max, startValue, AllowedSolution.ANY_SIDE);
    126        }
    127
    128        /**
    129         * {@inheritDoc}
    130         *
    131         * @throws ConvergenceException if the algorithm failed due to finite
    132         * precision.
    133         */
    134        @Override
    135        protected final double doSolve()
    136            throws ConvergenceException,
    137                   MathInternalError {
    138            // Get initial solution
    139            double x0 = getMin();
    140            double x1 = getMax();
    141            double f0 = computeObjectiveValue(x0);
    142            double f1 = computeObjectiveValue(x1);
    143
    144            // If one of the bounds is the exact root, return it. Since these are
    145            // not under-approximations or over-approximations, we can return them
    146            // regardless of the allowed solutions.
    147            if (f0 == 0.0) {
    148                return x0;
    149            }
    150            if (f1 == 0.0) {
    151                return x1;
    152            }
    153
    154            // Verify bracketing of initial solution.
    155            verifyBracketing(x0, x1);
    156
    157            // Get accuracies.
    158            final double ftol = getFunctionValueAccuracy();
    159            final double atol = getAbsoluteAccuracy();
    160            final double rtol = getRelativeAccuracy();
    161
    162            // Keep track of inverted intervals, meaning that the left bound is
    163            // larger than the right bound.
    164            boolean inverted = false;
    165
    166            // Keep finding better approximations.
    167            while (true) {
    168                // Calculate the next approximation.
    169                final double x = x1 - ((f1 * (x1 - x0)) / (f1 - f0));
    170                final double fx = computeObjectiveValue(x);
    171
    172                // If the new approximation is the exact root, return it. Since
    173                // this is not an under-approximation or an over-approximation,
    174                // we can return it regardless of the allowed solutions.
    175                if (fx == 0.0) {
    176                    return x;
    177                }
    178
    179                // Update the bounds with the new approximation.
    180                if (f1 * fx < 0) {
    181                    // The value of x1 has switched to the other bound, thus inverting
    182                    // the interval.
    183                    x0 = x1;
    184                    f0 = f1;
    185                    inverted = !inverted;
    186                } else {
    187                    switch (method) {
    188                    case ILLINOIS:
    189                        f0 *= 0.5;
    190                        break;
    191                    case PEGASUS:
    192                        f0 *= f1 / (f1 + fx);
    193                        break;
    194                    case REGULA_FALSI:
    195                        // Detect early that algorithm is stuck, instead of waiting
    196                        // for the maximum number of iterations to be exceeded.
    197                        if (x == x1) {
    198                            throw new ConvergenceException();
    199                        }
    200                        break;
    201                    default:
    202                        // Should never happen.
    203                        throw new MathInternalError();
    204                    }
    205                }
    206                // Update from [x0, x1] to [x0, x].
    207                x1 = x;
    208                f1 = fx;
    209
    210                // If the function value of the last approximation is too small,
    211                // given the function value accuracy, then we can't get closer to
    212                // the root than we already are.
    213                if (FastMath.abs(f1) <= ftol) {
    214                    switch (allowed) {
    215                    case ANY_SIDE:
    216                        return x1;
    217                    case LEFT_SIDE:
    218                        if (inverted) {
    219                            return x1;
    220                        }
    221                        break;
    222                    case RIGHT_SIDE:
    223                        if (!inverted) {
    224                            return x1;
    225                        }
    226                        break;
    227                    case BELOW_SIDE:
    228                        if (f1 <= 0) {
    229                            return x1;
    230                        }
    231                        break;
    232                    case ABOVE_SIDE:
    233                        if (f1 >= 0) {
    234                            return x1;
    235                        }
    236                        break;
    237                    default:
    238                        throw new MathInternalError();
    239                    }
    240                }
    241
    242                // If the current interval is within the given accuracies, we
    243                // are satisfied with the current approximation.
    244                if (FastMath.abs(x1 - x0) < FastMath.max(rtol * FastMath.abs(x1),
    245                                                         atol)) {
    246                    switch (allowed) {
    247                    case ANY_SIDE:
    248                        return x1;
    249                    case LEFT_SIDE:
    250                        return inverted ? x1 : x0;
    251                    case RIGHT_SIDE:
    252                        return inverted ? x0 : x1;
    253                    case BELOW_SIDE:
    254                        return (f1 <= 0) ? x1 : x0;
    255                    case ABOVE_SIDE:
    256                        return (f1 >= 0) ? x1 : x0;
    257                    default:
    258                        throw new MathInternalError();
    259                    }
    260                }
    261            }
    262        }
    263
    264        /** <em>Secant</em>-based root-finding methods. */
    265        protected enum Method {
    266
    267            /**
    268             * The {@link RegulaFalsiSolver <em>Regula Falsi</em>} or
    269             * <em>False Position</em> method.
    270             */
    271            REGULA_FALSI,
    272
    273            /** The {@link IllinoisSolver <em>Illinois</em>} method. */
    274            ILLINOIS,
    275
    276            /** The {@link PegasusSolver <em>Pegasus</em>} method. */
    277            PEGASUS;
    278
    279        }
    280    }