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.math4.legacy.analysis.solvers;
18
19 import org.apache.commons.math4.legacy.exception.NoBracketingException;
20 import org.apache.commons.math4.legacy.exception.TooManyEvaluationsException;
21 import org.apache.commons.math4.core.jdkmath.JdkMath;
22
23 /**
24 * Implements the <a href="http://mathworld.wolfram.com/RiddersMethod.html">
25 * Ridders' Method</a> for root finding of real univariate functions. For
26 * reference, see C. Ridders, <i>A new algorithm for computing a single root
27 * of a real continuous function </i>, IEEE Transactions on Circuits and
28 * Systems, 26 (1979), 979 - 980.
29 * <p>
30 * The function should be continuous but not necessarily smooth.</p>
31 *
32 * @since 1.2
33 */
34 public class RiddersSolver extends AbstractUnivariateSolver {
35 /** Default absolute accuracy. */
36 private static final double DEFAULT_ABSOLUTE_ACCURACY = 1e-6;
37
38 /**
39 * Construct a solver with default accuracy (1e-6).
40 */
41 public RiddersSolver() {
42 this(DEFAULT_ABSOLUTE_ACCURACY);
43 }
44 /**
45 * Construct a solver.
46 *
47 * @param absoluteAccuracy Absolute accuracy.
48 */
49 public RiddersSolver(double absoluteAccuracy) {
50 super(absoluteAccuracy);
51 }
52 /**
53 * Construct a solver.
54 *
55 * @param relativeAccuracy Relative accuracy.
56 * @param absoluteAccuracy Absolute accuracy.
57 */
58 public RiddersSolver(double relativeAccuracy,
59 double absoluteAccuracy) {
60 super(relativeAccuracy, absoluteAccuracy);
61 }
62
63 /**
64 * {@inheritDoc}
65 */
66 @Override
67 protected double doSolve()
68 throws TooManyEvaluationsException,
69 NoBracketingException {
70 double min = getMin();
71 double max = getMax();
72 // [x1, x2] is the bracketing interval in each iteration
73 // x3 is the midpoint of [x1, x2]
74 // x is the new root approximation and an endpoint of the new interval
75 double x1 = min;
76 double y1 = computeObjectiveValue(x1);
77 double x2 = max;
78 double y2 = computeObjectiveValue(x2);
79
80 // check for zeros before verifying bracketing
81 if (y1 == 0) {
82 return min;
83 }
84 if (y2 == 0) {
85 return max;
86 }
87 verifyBracketing(min, max);
88
89 final double absoluteAccuracy = getAbsoluteAccuracy();
90 final double functionValueAccuracy = getFunctionValueAccuracy();
91 final double relativeAccuracy = getRelativeAccuracy();
92
93 double oldx = Double.POSITIVE_INFINITY;
94 while (true) {
95 // calculate the new root approximation
96 final double x3 = 0.5 * (x1 + x2);
97 final double y3 = computeObjectiveValue(x3);
98 if (JdkMath.abs(y3) <= functionValueAccuracy) {
99 return x3;
100 }
101 final double delta = 1 - (y1 * y2) / (y3 * y3); // delta > 1 due to bracketing
102 final double correction = (JdkMath.signum(y2) * JdkMath.signum(y3)) *
103 (x3 - x1) / JdkMath.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 = JdkMath.max(relativeAccuracy * JdkMath.abs(x), absoluteAccuracy);
109 if (JdkMath.abs(x - oldx) <= tolerance) {
110 return x;
111 }
112 if (JdkMath.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 (JdkMath.signum(y1) + JdkMath.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 (JdkMath.signum(y2) + JdkMath.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 }