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 * https://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
18 /*
19 * This is not the original file distributed by the Apache Software Foundation
20 * It has been modified by the Hipparchus project
21 */
22 package org.hipparchus.analysis.solvers;
23
24 import org.hipparchus.exception.MathIllegalArgumentException;
25 import org.hipparchus.exception.MathIllegalStateException;
26 import org.hipparchus.util.FastMath;
27
28 /**
29 * Implements the <a href="http://mathworld.wolfram.com/RiddersMethod.html">
30 * Ridders' Method</a> for root finding of real univariate functions. For
31 * reference, see C. Ridders, <i>A new algorithm for computing a single root
32 * of a real continuous function </i>, IEEE Transactions on Circuits and
33 * Systems, 26 (1979), 979 - 980.
34 * <p>
35 * The function should be continuous but not necessarily smooth.</p>
36 *
37 */
38 public class RiddersSolver extends AbstractUnivariateSolver {
39 /** Default absolute accuracy. */
40 private static final double DEFAULT_ABSOLUTE_ACCURACY = 1e-6;
41
42 /**
43 * Construct a solver with default accuracy (1e-6).
44 */
45 public RiddersSolver() {
46 this(DEFAULT_ABSOLUTE_ACCURACY);
47 }
48 /**
49 * Construct a solver.
50 *
51 * @param absoluteAccuracy Absolute accuracy.
52 */
53 public RiddersSolver(double absoluteAccuracy) {
54 super(absoluteAccuracy);
55 }
56 /**
57 * Construct a solver.
58 *
59 * @param relativeAccuracy Relative accuracy.
60 * @param absoluteAccuracy Absolute accuracy.
61 */
62 public RiddersSolver(double relativeAccuracy,
63 double absoluteAccuracy) {
64 super(relativeAccuracy, absoluteAccuracy);
65 }
66
67 /**
68 * {@inheritDoc}
69 */
70 @Override
71 protected double doSolve()
72 throws MathIllegalArgumentException, MathIllegalStateException {
73 double min = getMin();
74 double max = getMax();
75 // [x1, x2] is the bracketing interval in each iteration
76 // x3 is the midpoint of [x1, x2]
77 // x is the new root approximation and an endpoint of the new interval
78 double x1 = min;
79 double y1 = computeObjectiveValue(x1);
80 double x2 = max;
81 double y2 = computeObjectiveValue(x2);
82
83 // check for zeros before verifying bracketing
84 if (y1 == 0) {
85 return min;
86 }
87 if (y2 == 0) {
88 return max;
89 }
90 verifyBracketing(min, max);
91
92 final double absoluteAccuracy = getAbsoluteAccuracy();
93 final double functionValueAccuracy = getFunctionValueAccuracy();
94 final double relativeAccuracy = getRelativeAccuracy();
95
96 double oldx = Double.POSITIVE_INFINITY;
97 while (true) {
98 // calculate the new root approximation
99 final double x3 = 0.5 * (x1 + x2);
100 final double y3 = computeObjectiveValue(x3);
101 if (FastMath.abs(y3) <= functionValueAccuracy) {
102 return x3;
103 }
104 final double delta = 1 - (y1 * y2) / (y3 * y3); // delta > 1 due to bracketing
105 final double correction = (FastMath.signum(y2) * FastMath.signum(y3)) *
106 (x3 - x1) / FastMath.sqrt(delta);
107 final double x = x3 - correction; // correction != 0
108 final double y = computeObjectiveValue(x);
109
110 // check for convergence
111 final double tolerance = FastMath.max(relativeAccuracy * FastMath.abs(x), absoluteAccuracy);
112 if (FastMath.abs(x - oldx) <= tolerance) {
113 return x;
114 }
115 if (FastMath.abs(y) <= functionValueAccuracy) {
116 return x;
117 }
118
119 // prepare the new interval for next iteration
120 // Ridders' method guarantees x1 < x < x2
121 if (correction > 0.0) { // x1 < x < x3
122 if (FastMath.signum(y1) + FastMath.signum(y) == 0.0) {
123 x2 = x;
124 y2 = y;
125 } else {
126 x1 = x;
127 x2 = x3;
128 y1 = y;
129 y2 = y3;
130 }
131 } else { // x3 < x < x2
132 if (FastMath.signum(y2) + FastMath.signum(y) == 0.0) {
133 x1 = x;
134 y1 = y;
135 } else {
136 x1 = x3;
137 x2 = x;
138 y1 = y3;
139 y2 = y;
140 }
141 }
142 oldx = x;
143 }
144 }
145 }