1 /*
2 * Licensed to the Apache Software Foundation (ASF) under one or more
3 * contributor license agreements. See the NOTICE file distributed with
4 * this work for additional information regarding copyright ownership.
5 * The ASF licenses this file to You under the Apache License, Version 2.0
6 * (the "License"); you may not use this file except in compliance with
7 * the License. You may obtain a copy of the License at
8 *
9 * http://www.apache.org/licenses/LICENSE-2.0
10 *
11 * Unless required by applicable law or agreed to in writing, software
12 * distributed under the License is distributed on an "AS IS" BASIS,
13 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14 * See the License for the specific language governing permissions and
15 * limitations under the License.
16 */
17 package org.apache.commons.math.analysis.solvers;
18
19 import org.apache.commons.math.util.FastMath;
20
21 /**
22 * Implements the <a href="http://mathworld.wolfram.com/RiddersMethod.html">
23 * Ridders' Method</a> for root finding of real univariate functions. For
24 * reference, see C. Ridders, <i>A new algorithm for computing a single root
25 * of a real continuous function </i>, IEEE Transactions on Circuits and
26 * Systems, 26 (1979), 979 - 980.
27 * <p>
28 * The function should be continuous but not necessarily smooth.</p>
29 *
30 * @version $Id: RiddersSolver.java 1183138 2011-10-13 22:21:04Z erans $
31 * @since 1.2
32 */
33 public class RiddersSolver extends AbstractUnivariateRealSolver {
34 /** Default absolute accuracy. */
35 private static final double DEFAULT_ABSOLUTE_ACCURACY = 1e-6;
36
37 /**
38 * Construct a solver with default accuracy (1e-6).
39 */
40 public RiddersSolver() {
41 this(DEFAULT_ABSOLUTE_ACCURACY);
42 }
43 /**
44 * Construct a solver.
45 *
46 * @param absoluteAccuracy Absolute accuracy.
47 */
48 public RiddersSolver(double absoluteAccuracy) {
49 super(absoluteAccuracy);
50 }
51 /**
52 * Construct a solver.
53 *
54 * @param relativeAccuracy Relative accuracy.
55 * @param absoluteAccuracy Absolute accuracy.
56 */
57 public RiddersSolver(double relativeAccuracy,
58 double absoluteAccuracy) {
59 super(relativeAccuracy, absoluteAccuracy);
60 }
61
62 /**
63 * {@inheritDoc}
64 */
65 @Override
66 protected double doSolve() {
67 double min = getMin();
68 double max = getMax();
69 // [x1, x2] is the bracketing interval in each iteration
70 // x3 is the midpoint of [x1, x2]
71 // x is the new root approximation and an endpoint of the new interval
72 double x1 = min;
73 double y1 = computeObjectiveValue(x1);
74 double x2 = max;
75 double y2 = computeObjectiveValue(x2);
76
77 // check for zeros before verifying bracketing
78 if (y1 == 0) {
79 return min;
80 }
81 if (y2 == 0) {
82 return max;
83 }
84 verifyBracketing(min, max);
85
86 final double absoluteAccuracy = getAbsoluteAccuracy();
87 final double functionValueAccuracy = getFunctionValueAccuracy();
88 final double relativeAccuracy = getRelativeAccuracy();
89
90 double oldx = Double.POSITIVE_INFINITY;
91 while (true) {
92 // calculate the new root approximation
93 final double x3 = 0.5 * (x1 + x2);
94 final double y3 = computeObjectiveValue(x3);
95 if (FastMath.abs(y3) <= functionValueAccuracy) {
96 return x3;
97 }
98 final double delta = 1 - (y1 * y2) / (y3 * y3); // delta > 1 due to bracketing
99 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 }