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    }