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 * @version $Id: RiddersSolver.java 1379560 2012-08-31 19:40:30Z erans $
033 * @since 1.2
034 */
035public class RiddersSolver extends AbstractUnivariateSolver {
036    /** Default absolute accuracy. */
037    private static final double DEFAULT_ABSOLUTE_ACCURACY = 1e-6;
038
039    /**
040     * Construct a solver with default accuracy (1e-6).
041     */
042    public RiddersSolver() {
043        this(DEFAULT_ABSOLUTE_ACCURACY);
044    }
045    /**
046     * Construct a solver.
047     *
048     * @param absoluteAccuracy Absolute accuracy.
049     */
050    public RiddersSolver(double absoluteAccuracy) {
051        super(absoluteAccuracy);
052    }
053    /**
054     * Construct a solver.
055     *
056     * @param relativeAccuracy Relative accuracy.
057     * @param absoluteAccuracy Absolute accuracy.
058     */
059    public RiddersSolver(double relativeAccuracy,
060                         double absoluteAccuracy) {
061        super(relativeAccuracy, absoluteAccuracy);
062    }
063
064    /**
065     * {@inheritDoc}
066     */
067    @Override
068    protected double doSolve()
069        throws TooManyEvaluationsException,
070               NoBracketingException {
071        double min = getMin();
072        double max = getMax();
073        // [x1, x2] is the bracketing interval in each iteration
074        // x3 is the midpoint of [x1, x2]
075        // x is the new root approximation and an endpoint of the new interval
076        double x1 = min;
077        double y1 = computeObjectiveValue(x1);
078        double x2 = max;
079        double y2 = computeObjectiveValue(x2);
080
081        // check for zeros before verifying bracketing
082        if (y1 == 0) {
083            return min;
084        }
085        if (y2 == 0) {
086            return max;
087        }
088        verifyBracketing(min, max);
089
090        final double absoluteAccuracy = getAbsoluteAccuracy();
091        final double functionValueAccuracy = getFunctionValueAccuracy();
092        final double relativeAccuracy = getRelativeAccuracy();
093
094        double oldx = Double.POSITIVE_INFINITY;
095        while (true) {
096            // calculate the new root approximation
097            final double x3 = 0.5 * (x1 + x2);
098            final double y3 = computeObjectiveValue(x3);
099            if (FastMath.abs(y3) <= functionValueAccuracy) {
100                return x3;
101            }
102            final double delta = 1 - (y1 * y2) / (y3 * y3);  // delta > 1 due to bracketing
103            final double correction = (FastMath.signum(y2) * FastMath.signum(y3)) *
104                                      (x3 - x1) / FastMath.sqrt(delta);
105            final double x = x3 - correction;                // correction != 0
106            final double y = computeObjectiveValue(x);
107
108            // check for convergence
109            final double tolerance = FastMath.max(relativeAccuracy * FastMath.abs(x), absoluteAccuracy);
110            if (FastMath.abs(x - oldx) <= tolerance) {
111                return x;
112            }
113            if (FastMath.abs(y) <= functionValueAccuracy) {
114                return x;
115            }
116
117            // prepare the new interval for next iteration
118            // Ridders' method guarantees x1 < x < x2
119            if (correction > 0.0) {             // x1 < x < x3
120                if (FastMath.signum(y1) + FastMath.signum(y) == 0.0) {
121                    x2 = x;
122                    y2 = y;
123                } else {
124                    x1 = x;
125                    x2 = x3;
126                    y1 = y;
127                    y2 = y3;
128                }
129            } else {                            // x3 < x < x2
130                if (FastMath.signum(y2) + FastMath.signum(y) == 0.0) {
131                    x1 = x;
132                    y1 = y;
133                } else {
134                    x1 = x3;
135                    x2 = x;
136                    y1 = y3;
137                    y2 = y;
138                }
139            }
140            oldx = x;
141        }
142    }
143}