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