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.math3.analysis.solvers;
018
019 import org.apache.commons.math3.util.FastMath;
020 import org.apache.commons.math3.exception.NoBracketingException;
021 import 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 */
035 public 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 }