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.analysis.UnivariateFunction;
26 import org.hipparchus.exception.LocalizedCoreFormats;
27 import org.hipparchus.exception.MathIllegalArgumentException;
28 import org.hipparchus.exception.MathIllegalStateException;
29 import org.hipparchus.exception.MathRuntimeException;
30 import org.hipparchus.util.FastMath;
31
32 /**
33 * Base class for all bracketing <em>Secant</em>-based methods for root-finding
34 * (approximating a zero of a univariate real function).
35 *
36 * <p>Implementation of the {@link RegulaFalsiSolver <em>Regula Falsi</em>} and
37 * {@link IllinoisSolver <em>Illinois</em>} methods is based on the
38 * following article: M. Dowell and P. Jarratt,
39 * <em>A modified regula falsi method for computing the root of an
40 * equation</em>, BIT Numerical Mathematics, volume 11, number 2,
41 * pages 168-174, Springer, 1971.</p>
42 *
43 * <p>Implementation of the {@link PegasusSolver <em>Pegasus</em>} method is
44 * based on the following article: M. Dowell and P. Jarratt,
45 * <em>The "Pegasus" method for computing the root of an equation</em>,
46 * BIT Numerical Mathematics, volume 12, number 4, pages 503-508, Springer,
47 * 1972.</p>
48 *
49 * <p>The {@link SecantSolver <em>Secant</em>} method is <em>not</em> a
50 * bracketing method, so it is not implemented here. It has a separate
51 * implementation.</p>
52 *
53 */
54 public abstract class BaseSecantSolver
55 extends AbstractUnivariateSolver
56 implements BracketedUnivariateSolver<UnivariateFunction> {
57
58 /** Default absolute accuracy. */
59 protected static final double DEFAULT_ABSOLUTE_ACCURACY = 1e-6;
60
61 /** The kinds of solutions that the algorithm may accept. */
62 private AllowedSolution allowed;
63
64 /** The <em>Secant</em>-based root-finding method to use. */
65 private final Method method;
66
67 /**
68 * Construct a solver.
69 *
70 * @param absoluteAccuracy Absolute accuracy.
71 * @param method <em>Secant</em>-based root-finding method to use.
72 */
73 protected BaseSecantSolver(final double absoluteAccuracy, final Method method) {
74 super(absoluteAccuracy);
75 this.allowed = AllowedSolution.ANY_SIDE;
76 this.method = method;
77 }
78
79 /**
80 * Construct a solver.
81 *
82 * @param relativeAccuracy Relative accuracy.
83 * @param absoluteAccuracy Absolute accuracy.
84 * @param method <em>Secant</em>-based root-finding method to use.
85 */
86 protected BaseSecantSolver(final double relativeAccuracy,
87 final double absoluteAccuracy,
88 final Method method) {
89 super(relativeAccuracy, absoluteAccuracy);
90 this.allowed = AllowedSolution.ANY_SIDE;
91 this.method = method;
92 }
93
94 /**
95 * Construct a solver.
96 *
97 * @param relativeAccuracy Maximum relative error.
98 * @param absoluteAccuracy Maximum absolute error.
99 * @param functionValueAccuracy Maximum function value error.
100 * @param method <em>Secant</em>-based root-finding method to use
101 */
102 protected BaseSecantSolver(final double relativeAccuracy,
103 final double absoluteAccuracy,
104 final double functionValueAccuracy,
105 final Method method) {
106 super(relativeAccuracy, absoluteAccuracy, functionValueAccuracy);
107 this.allowed = AllowedSolution.ANY_SIDE;
108 this.method = method;
109 }
110
111 /** {@inheritDoc} */
112 @Override
113 public double solve(final int maxEval, final UnivariateFunction f,
114 final double min, final double max,
115 final AllowedSolution allowedSolution) {
116 return solve(maxEval, f, min, max, min + 0.5 * (max - min), allowedSolution);
117 }
118
119 /** {@inheritDoc} */
120 @Override
121 public double solve(final int maxEval, final UnivariateFunction f,
122 final double min, final double max, final double startValue,
123 final AllowedSolution allowedSolution) {
124 this.allowed = allowedSolution;
125 return super.solve(maxEval, f, min, max, startValue);
126 }
127
128 /** {@inheritDoc} */
129 @Override
130 public double solve(final int maxEval, final UnivariateFunction f,
131 final double min, final double max, final double startValue) {
132 return solve(maxEval, f, min, max, startValue, AllowedSolution.ANY_SIDE);
133 }
134
135 /** {@inheritDoc} */
136 @Override
137 public Interval solveInterval(final int maxEval,
138 final UnivariateFunction f,
139 final double min,
140 final double max,
141 final double startValue) throws MathIllegalArgumentException, MathIllegalStateException {
142 setup(maxEval, f, min, max, startValue);
143 this.allowed = null;
144 return doSolveInterval();
145 }
146
147 /**
148 * {@inheritDoc}
149 *
150 * @throws MathIllegalStateException if the algorithm failed due to finite
151 * precision.
152 */
153 @Override
154 protected final double doSolve() throws MathIllegalStateException {
155 return doSolveInterval().getSide(allowed);
156 }
157
158 /**
159 * Find a root and return the containing interval.
160 *
161 * @return an interval containing the root such that the selected end point meets the
162 * convergence criteria.
163 * @throws MathIllegalStateException if convergence fails.
164 */
165 protected final Interval doSolveInterval() throws MathIllegalStateException {
166 // Get initial solution
167 double x0 = getMin();
168 double x1 = getMax();
169 double f0 = computeObjectiveValue(x0);
170 double f1 = computeObjectiveValue(x1);
171
172 // If one of the bounds is the exact root, return it. Since these are
173 // not under-approximations or over-approximations, we can return them
174 // regardless of the allowed solutions.
175 if (f0 == 0.0) {
176 return new Interval(x0, f0, x0, f0);
177 }
178 if (f1 == 0.0) {
179 return new Interval(x1, f1, x1, f1);
180 }
181
182 // Verify bracketing of initial solution.
183 verifyBracketing(x0, x1);
184
185 // Get accuracies.
186 final double ftol = getFunctionValueAccuracy();
187 final double atol = getAbsoluteAccuracy();
188 final double rtol = getRelativeAccuracy();
189
190 // Keep track of inverted intervals, meaning that the left bound is
191 // larger than the right bound.
192 boolean inverted = false;
193
194 // Keep finding better approximations.
195 while (true) {
196 // Calculate the next approximation.
197 final double x = x1 - ((f1 * (x1 - x0)) / (f1 - f0));
198 final double fx = computeObjectiveValue(x);
199
200 // If the new approximation is the exact root, return it. Since
201 // this is not an under-approximation or an over-approximation,
202 // we can return it regardless of the allowed solutions.
203 if (fx == 0.0) {
204 return new Interval(x, fx, x, fx);
205 }
206
207 // Update the bounds with the new approximation.
208 if (f1 * fx < 0) {
209 // The value of x1 has switched to the other bound, thus inverting
210 // the interval.
211 x0 = x1;
212 f0 = f1;
213 inverted = !inverted;
214 } else {
215 switch (method) {
216 case ILLINOIS:
217 f0 *= 0.5;
218 break;
219 case PEGASUS:
220 f0 *= f1 / (f1 + fx);
221 break;
222 case REGULA_FALSI:
223 // Detect early that algorithm is stuck, instead of waiting
224 // for the maximum number of iterations to be exceeded.
225 if (x == x1) {
226 throw new MathIllegalStateException(LocalizedCoreFormats.CONVERGENCE_FAILED);
227 }
228 break;
229 default:
230 // Should never happen.
231 throw MathRuntimeException.createInternalError();
232 }
233 }
234 // Update from [x0, x1] to [x0, x].
235 x1 = x;
236 f1 = fx;
237
238 // If the current interval is within the given accuracies, we
239 // are satisfied with the current approximation.
240 if (FastMath.abs(x1 - x0) < FastMath.max(rtol * FastMath.abs(x1), atol) ||
241 (FastMath.abs(f1) < ftol && (allowed == AllowedSolution.ANY_SIDE ||
242 (inverted && allowed == AllowedSolution.LEFT_SIDE) ||
243 (!inverted && allowed == AllowedSolution.RIGHT_SIDE) ||
244 (f1 <= 0.0 && allowed == AllowedSolution.BELOW_SIDE) ||
245 (f1 >= 0.0 && allowed == AllowedSolution.ABOVE_SIDE)))) {
246 if (inverted) {
247 return new Interval(x1, f1, x0, f0);
248 } else {
249 return new Interval(x0, f0, x1, f1);
250 }
251 }
252 }
253 }
254
255 /** <em>Secant</em>-based root-finding methods. */
256 protected enum Method {
257
258 /**
259 * The {@link RegulaFalsiSolver <em>Regula Falsi</em>} or
260 * <em>False Position</em> method.
261 */
262 REGULA_FALSI,
263
264 /** The {@link IllinoisSolver <em>Illinois</em>} method. */
265 ILLINOIS,
266
267 /** The {@link PegasusSolver <em>Pegasus</em>} method. */
268 PEGASUS;
269
270 }
271 }