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
23 package org.hipparchus.analysis.solvers;
24
25 import org.hipparchus.exception.MathIllegalArgumentException;
26 import org.hipparchus.exception.MathIllegalStateException;
27 import org.hipparchus.util.FastMath;
28
29 /**
30 * Implements the <em>Secant</em> method for root-finding (approximating a
31 * zero of a univariate real function). The solution that is maintained is
32 * not bracketed, and as such convergence is not guaranteed.
33 *
34 * <p>Implementation based on the following article: M. Dowell and P. Jarratt,
35 * <em>A modified regula falsi method for computing the root of an
36 * equation</em>, BIT Numerical Mathematics, volume 11, number 2,
37 * pages 168-174, Springer, 1971.</p>
38 *
39 * <p>Note that since release 3.0 this class implements the actual
40 * <em>Secant</em> algorithm, and not a modified one. As such, the 3.0 version
41 * is not backwards compatible with previous versions. To use an algorithm
42 * similar to the pre-3.0 releases, use the
43 * {@link IllinoisSolver <em>Illinois</em>} algorithm or the
44 * {@link PegasusSolver <em>Pegasus</em>} algorithm.</p>
45 *
46 */
47 public class SecantSolver extends AbstractUnivariateSolver {
48
49 /** Default absolute accuracy. */
50 protected static final double DEFAULT_ABSOLUTE_ACCURACY = 1e-6;
51
52 /** Construct a solver with default accuracy (1e-6). */
53 public SecantSolver() {
54 super(DEFAULT_ABSOLUTE_ACCURACY);
55 }
56
57 /**
58 * Construct a solver.
59 *
60 * @param absoluteAccuracy absolute accuracy
61 */
62 public SecantSolver(final double absoluteAccuracy) {
63 super(absoluteAccuracy);
64 }
65
66 /**
67 * Construct a solver.
68 *
69 * @param relativeAccuracy relative accuracy
70 * @param absoluteAccuracy absolute accuracy
71 */
72 public SecantSolver(final double relativeAccuracy,
73 final double absoluteAccuracy) {
74 super(relativeAccuracy, absoluteAccuracy);
75 }
76
77 /** {@inheritDoc} */
78 @Override
79 protected final double doSolve()
80 throws MathIllegalArgumentException, MathIllegalStateException {
81 // Get initial solution
82 double x0 = getMin();
83 double x1 = getMax();
84 double f0 = computeObjectiveValue(x0);
85 double f1 = computeObjectiveValue(x1);
86
87 // If one of the bounds is the exact root, return it. Since these are
88 // not under-approximations or over-approximations, we can return them
89 // regardless of the allowed solutions.
90 if (f0 == 0.0) {
91 return x0;
92 }
93 if (f1 == 0.0) {
94 return x1;
95 }
96
97 // Verify bracketing of initial solution.
98 verifyBracketing(x0, x1);
99
100 // Get accuracies.
101 final double ftol = getFunctionValueAccuracy();
102 final double atol = getAbsoluteAccuracy();
103 final double rtol = getRelativeAccuracy();
104
105 // Keep finding better approximations.
106 while (true) {
107 // Calculate the next approximation.
108 final double x = x1 - ((f1 * (x1 - x0)) / (f1 - f0));
109 final double fx = computeObjectiveValue(x);
110
111 // If the new approximation is the exact root, return it. Since
112 // this is not an under-approximation or an over-approximation,
113 // we can return it regardless of the allowed solutions.
114 if (fx == 0.0) {
115 return x;
116 }
117
118 // Update the bounds with the new approximation.
119 x0 = x1;
120 f0 = f1;
121 x1 = x;
122 f1 = fx;
123
124 // If the function value of the last approximation is too small,
125 // given the function value accuracy, then we can't get closer to
126 // the root than we already are.
127 if (FastMath.abs(f1) <= ftol) {
128 return x1;
129 }
130
131 // If the current interval is within the given accuracies, we
132 // are satisfied with the current approximation.
133 if (FastMath.abs(x1 - x0) < FastMath.max(rtol * FastMath.abs(x1), atol)) {
134 return x1;
135 }
136 }
137 }
138
139 }