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