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.math4.analysis.solvers;
019
020import org.apache.commons.math4.analysis.UnivariateFunction;
021import org.apache.commons.math4.exception.ConvergenceException;
022import org.apache.commons.math4.exception.MathInternalError;
023import org.apache.commons.math4.util.FastMath;
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    @Override
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    @Override
115    public double solve(final int maxEval, final UnivariateFunction f,
116                        final double min, final double max, final double startValue,
117                        final AllowedSolution allowedSolution) {
118        this.allowed = allowedSolution;
119        return super.solve(maxEval, f, min, max, startValue);
120    }
121
122    /** {@inheritDoc} */
123    @Override
124    public double solve(final int maxEval, final UnivariateFunction f,
125                        final double min, final double max, final double startValue) {
126        return solve(maxEval, f, min, max, startValue, AllowedSolution.ANY_SIDE);
127    }
128
129    /**
130     * {@inheritDoc}
131     *
132     * @throws ConvergenceException if the algorithm failed due to finite
133     * precision.
134     */
135    @Override
136    protected final double doSolve()
137        throws ConvergenceException {
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}