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    package org.apache.commons.math.analysis.solvers;
018    
019    import org.apache.commons.math.ConvergenceException;
020    import org.apache.commons.math.FunctionEvaluationException;
021    import org.apache.commons.math.MaxIterationsExceededException;
022    import org.apache.commons.math.analysis.UnivariateRealFunction;
023    import org.apache.commons.math.util.FastMath;
024    import org.apache.commons.math.util.MathUtils;
025    
026    /**
027     * Implements the <a href="http://mathworld.wolfram.com/RiddersMethod.html">
028     * Ridders' Method</a> for root finding of real univariate functions. For
029     * reference, see C. Ridders, <i>A new algorithm for computing a single root
030     * of a real continuous function </i>, IEEE Transactions on Circuits and
031     * Systems, 26 (1979), 979 - 980.
032     * <p>
033     * The function should be continuous but not necessarily smooth.</p>
034     *
035     * @version $Revision: 1070725 $ $Date: 2011-02-15 02:31:12 +0100 (mar. 15 f??vr. 2011) $
036     * @since 1.2
037     */
038    public class RiddersSolver extends UnivariateRealSolverImpl {
039    
040        /**
041         * Construct a solver for the given function.
042         *
043         * @param f function to solve
044         * @deprecated as of 2.0 the function to solve is passed as an argument
045         * to the {@link #solve(UnivariateRealFunction, double, double)} or
046         * {@link UnivariateRealSolverImpl#solve(UnivariateRealFunction, double, double, double)}
047         * method.
048         */
049        @Deprecated
050        public RiddersSolver(UnivariateRealFunction f) {
051            super(f, 100, 1E-6);
052        }
053    
054        /**
055         * Construct a solver.
056         * @deprecated in 2.2
057         */
058        @Deprecated
059        public RiddersSolver() {
060            super(100, 1E-6);
061        }
062    
063        /** {@inheritDoc} */
064        @Deprecated
065        public double solve(final double min, final double max)
066            throws ConvergenceException, FunctionEvaluationException {
067            return solve(f, min, max);
068        }
069    
070        /** {@inheritDoc} */
071        @Deprecated
072        public double solve(final double min, final double max, final double initial)
073            throws ConvergenceException, FunctionEvaluationException {
074            return solve(f, min, max, initial);
075        }
076    
077        /**
078         * Find a root in the given interval with initial value.
079         * <p>
080         * Requires bracketing condition.</p>
081         *
082         * @param f the function to solve
083         * @param min the lower bound for the interval
084         * @param max the upper bound for the interval
085         * @param initial the start value to use
086         * @param maxEval Maximum number of evaluations.
087         * @return the point at which the function value is zero
088         * @throws MaxIterationsExceededException if the maximum iteration count is exceeded
089         * @throws FunctionEvaluationException if an error occurs evaluating the function
090         * @throws IllegalArgumentException if any parameters are invalid
091         */
092        @Override
093        public double solve(int maxEval, final UnivariateRealFunction f,
094                            final double min, final double max, final double initial)
095            throws MaxIterationsExceededException, FunctionEvaluationException {
096            setMaximalIterationCount(maxEval);
097            return solve(f, min, max, initial);
098        }
099    
100        /**
101         * Find a root in the given interval with initial value.
102         * <p>
103         * Requires bracketing condition.</p>
104         *
105         * @param f the function to solve
106         * @param min the lower bound for the interval
107         * @param max the upper bound for the interval
108         * @param initial the start value to use
109         * @return the point at which the function value is zero
110         * @throws MaxIterationsExceededException if the maximum iteration count is exceeded
111         * @throws FunctionEvaluationException if an error occurs evaluating the function
112         * @throws IllegalArgumentException if any parameters are invalid
113         * @deprecated in 2.2 (to be removed in 3.0).
114         */
115        @Deprecated
116        public double solve(final UnivariateRealFunction f,
117                            final double min, final double max, final double initial)
118            throws MaxIterationsExceededException, FunctionEvaluationException {
119    
120            // check for zeros before verifying bracketing
121            if (f.value(min) == 0.0) { return min; }
122            if (f.value(max) == 0.0) { return max; }
123            if (f.value(initial) == 0.0) { return initial; }
124    
125            verifyBracketing(min, max, f);
126            verifySequence(min, initial, max);
127            if (isBracketing(min, initial, f)) {
128                return solve(f, min, initial);
129            } else {
130                return solve(f, initial, max);
131            }
132        }
133    
134        /**
135         * Find a root in the given interval.
136         * <p>
137         * Requires bracketing condition.</p>
138         *
139         * @param f the function to solve
140         * @param min the lower bound for the interval
141         * @param max the upper bound for the interval
142         * @param maxEval Maximum number of evaluations.
143         * @return the point at which the function value is zero
144         * @throws MaxIterationsExceededException if the maximum iteration count is exceeded
145         * @throws FunctionEvaluationException if an error occurs evaluating the function
146         * @throws IllegalArgumentException if any parameters are invalid
147         */
148        @Override
149        public double solve(int maxEval, final UnivariateRealFunction f,
150                            final double min, final double max)
151            throws MaxIterationsExceededException, FunctionEvaluationException {
152            setMaximalIterationCount(maxEval);
153            return solve(f, min, max);
154        }
155    
156        /**
157         * Find a root in the given interval.
158         * <p>
159         * Requires bracketing condition.</p>
160         *
161         * @param f the function to solve
162         * @param min the lower bound for the interval
163         * @param max the upper bound for the interval
164         * @return the point at which the function value is zero
165         * @throws MaxIterationsExceededException if the maximum iteration count is exceeded
166         * @throws FunctionEvaluationException if an error occurs evaluating the function
167         * @throws IllegalArgumentException if any parameters are invalid
168         * @deprecated in 2.2 (to be removed in 3.0).
169         */
170        @Deprecated
171        public double solve(final UnivariateRealFunction f,
172                            final double min, final double max)
173            throws MaxIterationsExceededException, FunctionEvaluationException {
174    
175            // [x1, x2] is the bracketing interval in each iteration
176            // x3 is the midpoint of [x1, x2]
177            // x is the new root approximation and an endpoint of the new interval
178            double x1 = min;
179            double y1 = f.value(x1);
180            double x2 = max;
181            double y2 = f.value(x2);
182    
183            // check for zeros before verifying bracketing
184            if (y1 == 0.0) {
185                return min;
186            }
187            if (y2 == 0.0) {
188                return max;
189            }
190            verifyBracketing(min, max, f);
191    
192            int i = 1;
193            double oldx = Double.POSITIVE_INFINITY;
194            while (i <= maximalIterationCount) {
195                // calculate the new root approximation
196                final double x3 = 0.5 * (x1 + x2);
197                final double y3 = f.value(x3);
198                if (FastMath.abs(y3) <= functionValueAccuracy) {
199                    setResult(x3, i);
200                    return result;
201                }
202                final double delta = 1 - (y1 * y2) / (y3 * y3);  // delta > 1 due to bracketing
203                final double correction = (MathUtils.sign(y2) * MathUtils.sign(y3)) *
204                                          (x3 - x1) / FastMath.sqrt(delta);
205                final double x = x3 - correction;                // correction != 0
206                final double y = f.value(x);
207    
208                // check for convergence
209                final double tolerance = FastMath.max(relativeAccuracy * FastMath.abs(x), absoluteAccuracy);
210                if (FastMath.abs(x - oldx) <= tolerance) {
211                    setResult(x, i);
212                    return result;
213                }
214                if (FastMath.abs(y) <= functionValueAccuracy) {
215                    setResult(x, i);
216                    return result;
217                }
218    
219                // prepare the new interval for next iteration
220                // Ridders' method guarantees x1 < x < x2
221                if (correction > 0.0) {             // x1 < x < x3
222                    if (MathUtils.sign(y1) + MathUtils.sign(y) == 0.0) {
223                        x2 = x;
224                        y2 = y;
225                    } else {
226                        x1 = x;
227                        x2 = x3;
228                        y1 = y;
229                        y2 = y3;
230                    }
231                } else {                            // x3 < x < x2
232                    if (MathUtils.sign(y2) + MathUtils.sign(y) == 0.0) {
233                        x1 = x;
234                        y1 = y;
235                    } else {
236                        x1 = x3;
237                        x2 = x;
238                        y1 = y3;
239                        y2 = y;
240                    }
241                }
242                oldx = x;
243                i++;
244            }
245            throw new MaxIterationsExceededException(maximalIterationCount);
246        }
247    }