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 */
017package org.apache.commons.math3.analysis.solvers;
018
019import org.apache.commons.math3.util.FastMath;
020import org.apache.commons.math3.exception.NoBracketingException;
021import org.apache.commons.math3.exception.TooManyEvaluationsException;
022
023/**
024 * Implements the <a href="http://mathworld.wolfram.com/RiddersMethod.html">
025 * Ridders' Method</a> for root finding of real univariate functions. For
026 * reference, see C. Ridders, <i>A new algorithm for computing a single root
027 * of a real continuous function </i>, IEEE Transactions on Circuits and
028 * Systems, 26 (1979), 979 - 980.
029 * <p>
030 * The function should be continuous but not necessarily smooth.</p>
031 *
032 * @since 1.2
033 */
034public class RiddersSolver extends AbstractUnivariateSolver {
035    /** Default absolute accuracy. */
036    private static final double DEFAULT_ABSOLUTE_ACCURACY = 1e-6;
037
038    /**
039     * Construct a solver with default accuracy (1e-6).
040     */
041    public RiddersSolver() {
042        this(DEFAULT_ABSOLUTE_ACCURACY);
043    }
044    /**
045     * Construct a solver.
046     *
047     * @param absoluteAccuracy Absolute accuracy.
048     */
049    public RiddersSolver(double absoluteAccuracy) {
050        super(absoluteAccuracy);
051    }
052    /**
053     * Construct a solver.
054     *
055     * @param relativeAccuracy Relative accuracy.
056     * @param absoluteAccuracy Absolute accuracy.
057     */
058    public RiddersSolver(double relativeAccuracy,
059                         double absoluteAccuracy) {
060        super(relativeAccuracy, absoluteAccuracy);
061    }
062
063    /**
064     * {@inheritDoc}
065     */
066    @Override
067    protected double doSolve()
068        throws TooManyEvaluationsException,
069               NoBracketingException {
070        double min = getMin();
071        double max = getMax();
072        // [x1, x2] is the bracketing interval in each iteration
073        // x3 is the midpoint of [x1, x2]
074        // x is the new root approximation and an endpoint of the new interval
075        double x1 = min;
076        double y1 = computeObjectiveValue(x1);
077        double x2 = max;
078        double y2 = computeObjectiveValue(x2);
079
080        // check for zeros before verifying bracketing
081        if (y1 == 0) {
082            return min;
083        }
084        if (y2 == 0) {
085            return max;
086        }
087        verifyBracketing(min, max);
088
089        final double absoluteAccuracy = getAbsoluteAccuracy();
090        final double functionValueAccuracy = getFunctionValueAccuracy();
091        final double relativeAccuracy = getRelativeAccuracy();
092
093        double oldx = Double.POSITIVE_INFINITY;
094        while (true) {
095            // calculate the new root approximation
096            final double x3 = 0.5 * (x1 + x2);
097            final double y3 = computeObjectiveValue(x3);
098            if (FastMath.abs(y3) <= functionValueAccuracy) {
099                return x3;
100            }
101            final double delta = 1 - (y1 * y2) / (y3 * y3);  // delta > 1 due to bracketing
102            final double correction = (FastMath.signum(y2) * FastMath.signum(y3)) *
103                                      (x3 - x1) / FastMath.sqrt(delta);
104            final double x = x3 - correction;                // correction != 0
105            final double y = computeObjectiveValue(x);
106
107            // check for convergence
108            final double tolerance = FastMath.max(relativeAccuracy * FastMath.abs(x), absoluteAccuracy);
109            if (FastMath.abs(x - oldx) <= tolerance) {
110                return x;
111            }
112            if (FastMath.abs(y) <= functionValueAccuracy) {
113                return x;
114            }
115
116            // prepare the new interval for next iteration
117            // Ridders' method guarantees x1 < x < x2
118            if (correction > 0.0) {             // x1 < x < x3
119                if (FastMath.signum(y1) + FastMath.signum(y) == 0.0) {
120                    x2 = x;
121                    y2 = y;
122                } else {
123                    x1 = x;
124                    x2 = x3;
125                    y1 = y;
126                    y2 = y3;
127                }
128            } else {                            // x3 < x < x2
129                if (FastMath.signum(y2) + FastMath.signum(y) == 0.0) {
130                    x1 = x;
131                    y1 = y;
132                } else {
133                    x1 = x3;
134                    x2 = x;
135                    y1 = y3;
136                    y2 = y;
137                }
138            }
139            oldx = x;
140        }
141    }
142}