001 /*
002 * Licensed to the Apache Software Foundation (ASF) under one or more
003 * contributor license agreements. See the NOTICE file distributed with
004 * this work for additional information regarding copyright ownership.
005 * The ASF licenses this file to You under the Apache License, Version 2.0
006 * (the "License"); you may not use this file except in compliance with
007 * the License. You may obtain a copy of the License at
008 *
009 * http://www.apache.org/licenses/LICENSE-2.0
010 *
011 * Unless required by applicable law or agreed to in writing, software
012 * distributed under the License is distributed on an "AS IS" BASIS,
013 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
014 * See the License for the specific language governing permissions and
015 * limitations under the License.
016 */
017
018 package org.apache.commons.math3.analysis.solvers;
019
020 import org.apache.commons.math3.util.FastMath;
021 import org.apache.commons.math3.analysis.UnivariateFunction;
022 import org.apache.commons.math3.exception.ConvergenceException;
023 import org.apache.commons.math3.exception.MathInternalError;
024
025 /**
026 * Base class for all bracketing <em>Secant</em>-based methods for root-finding
027 * (approximating a zero of a univariate real function).
028 *
029 * <p>Implementation of the {@link RegulaFalsiSolver <em>Regula Falsi</em>} and
030 * {@link IllinoisSolver <em>Illinois</em>} methods is based on the
031 * following article: M. Dowell and P. Jarratt,
032 * <em>A modified regula falsi method for computing the root of an
033 * equation</em>, BIT Numerical Mathematics, volume 11, number 2,
034 * pages 168-174, Springer, 1971.</p>
035 *
036 * <p>Implementation of the {@link PegasusSolver <em>Pegasus</em>} method is
037 * based on the following article: M. Dowell and P. Jarratt,
038 * <em>The "Pegasus" method for computing the root of an equation</em>,
039 * BIT Numerical Mathematics, volume 12, number 4, pages 503-508, Springer,
040 * 1972.</p>
041 *
042 * <p>The {@link SecantSolver <em>Secant</em>} method is <em>not</em> a
043 * bracketing method, so it is not implemented here. It has a separate
044 * implementation.</p>
045 *
046 * @since 3.0
047 * @version $Id: BaseSecantSolver.java 1379560 2012-08-31 19:40:30Z erans $
048 */
049 public abstract class BaseSecantSolver
050 extends AbstractUnivariateSolver
051 implements BracketedUnivariateSolver<UnivariateFunction> {
052
053 /** Default absolute accuracy. */
054 protected static final double DEFAULT_ABSOLUTE_ACCURACY = 1e-6;
055
056 /** The kinds of solutions that the algorithm may accept. */
057 private AllowedSolution allowed;
058
059 /** The <em>Secant</em>-based root-finding method to use. */
060 private final Method method;
061
062 /**
063 * Construct a solver.
064 *
065 * @param absoluteAccuracy Absolute accuracy.
066 * @param method <em>Secant</em>-based root-finding method to use.
067 */
068 protected BaseSecantSolver(final double absoluteAccuracy, final Method method) {
069 super(absoluteAccuracy);
070 this.allowed = AllowedSolution.ANY_SIDE;
071 this.method = method;
072 }
073
074 /**
075 * Construct a solver.
076 *
077 * @param relativeAccuracy Relative accuracy.
078 * @param absoluteAccuracy Absolute accuracy.
079 * @param method <em>Secant</em>-based root-finding method to use.
080 */
081 protected BaseSecantSolver(final double relativeAccuracy,
082 final double absoluteAccuracy,
083 final Method method) {
084 super(relativeAccuracy, absoluteAccuracy);
085 this.allowed = AllowedSolution.ANY_SIDE;
086 this.method = method;
087 }
088
089 /**
090 * Construct a solver.
091 *
092 * @param relativeAccuracy Maximum relative error.
093 * @param absoluteAccuracy Maximum absolute error.
094 * @param functionValueAccuracy Maximum function value error.
095 * @param method <em>Secant</em>-based root-finding method to use
096 */
097 protected BaseSecantSolver(final double relativeAccuracy,
098 final double absoluteAccuracy,
099 final double functionValueAccuracy,
100 final Method method) {
101 super(relativeAccuracy, absoluteAccuracy, functionValueAccuracy);
102 this.allowed = AllowedSolution.ANY_SIDE;
103 this.method = method;
104 }
105
106 /** {@inheritDoc} */
107 public double solve(final int maxEval, final UnivariateFunction f,
108 final double min, final double max,
109 final AllowedSolution allowedSolution) {
110 return solve(maxEval, f, min, max, min + 0.5 * (max - min), allowedSolution);
111 }
112
113 /** {@inheritDoc} */
114 public double solve(final int maxEval, final UnivariateFunction f,
115 final double min, final double max, final double startValue,
116 final AllowedSolution allowedSolution) {
117 this.allowed = allowedSolution;
118 return super.solve(maxEval, f, min, max, startValue);
119 }
120
121 /** {@inheritDoc} */
122 @Override
123 public double solve(final int maxEval, final UnivariateFunction f,
124 final double min, final double max, final double startValue) {
125 return solve(maxEval, f, min, max, startValue, AllowedSolution.ANY_SIDE);
126 }
127
128 /**
129 * {@inheritDoc}
130 *
131 * @throws ConvergenceException if the algorithm failed due to finite
132 * precision.
133 */
134 @Override
135 protected final double doSolve()
136 throws ConvergenceException,
137 MathInternalError {
138 // Get initial solution
139 double x0 = getMin();
140 double x1 = getMax();
141 double f0 = computeObjectiveValue(x0);
142 double f1 = computeObjectiveValue(x1);
143
144 // If one of the bounds is the exact root, return it. Since these are
145 // not under-approximations or over-approximations, we can return them
146 // regardless of the allowed solutions.
147 if (f0 == 0.0) {
148 return x0;
149 }
150 if (f1 == 0.0) {
151 return x1;
152 }
153
154 // Verify bracketing of initial solution.
155 verifyBracketing(x0, x1);
156
157 // Get accuracies.
158 final double ftol = getFunctionValueAccuracy();
159 final double atol = getAbsoluteAccuracy();
160 final double rtol = getRelativeAccuracy();
161
162 // Keep track of inverted intervals, meaning that the left bound is
163 // larger than the right bound.
164 boolean inverted = false;
165
166 // Keep finding better approximations.
167 while (true) {
168 // Calculate the next approximation.
169 final double x = x1 - ((f1 * (x1 - x0)) / (f1 - f0));
170 final double fx = computeObjectiveValue(x);
171
172 // If the new approximation is the exact root, return it. Since
173 // this is not an under-approximation or an over-approximation,
174 // we can return it regardless of the allowed solutions.
175 if (fx == 0.0) {
176 return x;
177 }
178
179 // Update the bounds with the new approximation.
180 if (f1 * fx < 0) {
181 // The value of x1 has switched to the other bound, thus inverting
182 // the interval.
183 x0 = x1;
184 f0 = f1;
185 inverted = !inverted;
186 } else {
187 switch (method) {
188 case ILLINOIS:
189 f0 *= 0.5;
190 break;
191 case PEGASUS:
192 f0 *= f1 / (f1 + fx);
193 break;
194 case REGULA_FALSI:
195 // Detect early that algorithm is stuck, instead of waiting
196 // for the maximum number of iterations to be exceeded.
197 if (x == x1) {
198 throw new ConvergenceException();
199 }
200 break;
201 default:
202 // Should never happen.
203 throw new MathInternalError();
204 }
205 }
206 // Update from [x0, x1] to [x0, x].
207 x1 = x;
208 f1 = fx;
209
210 // If the function value of the last approximation is too small,
211 // given the function value accuracy, then we can't get closer to
212 // the root than we already are.
213 if (FastMath.abs(f1) <= ftol) {
214 switch (allowed) {
215 case ANY_SIDE:
216 return x1;
217 case LEFT_SIDE:
218 if (inverted) {
219 return x1;
220 }
221 break;
222 case RIGHT_SIDE:
223 if (!inverted) {
224 return x1;
225 }
226 break;
227 case BELOW_SIDE:
228 if (f1 <= 0) {
229 return x1;
230 }
231 break;
232 case ABOVE_SIDE:
233 if (f1 >= 0) {
234 return x1;
235 }
236 break;
237 default:
238 throw new MathInternalError();
239 }
240 }
241
242 // If the current interval is within the given accuracies, we
243 // are satisfied with the current approximation.
244 if (FastMath.abs(x1 - x0) < FastMath.max(rtol * FastMath.abs(x1),
245 atol)) {
246 switch (allowed) {
247 case ANY_SIDE:
248 return x1;
249 case LEFT_SIDE:
250 return inverted ? x1 : x0;
251 case RIGHT_SIDE:
252 return inverted ? x0 : x1;
253 case BELOW_SIDE:
254 return (f1 <= 0) ? x1 : x0;
255 case ABOVE_SIDE:
256 return (f1 >= 0) ? x1 : x0;
257 default:
258 throw new MathInternalError();
259 }
260 }
261 }
262 }
263
264 /** <em>Secant</em>-based root-finding methods. */
265 protected enum Method {
266
267 /**
268 * The {@link RegulaFalsiSolver <em>Regula Falsi</em>} or
269 * <em>False Position</em> method.
270 */
271 REGULA_FALSI,
272
273 /** The {@link IllinoisSolver <em>Illinois</em>} method. */
274 ILLINOIS,
275
276 /** The {@link PegasusSolver <em>Pegasus</em>} method. */
277 PEGASUS;
278
279 }
280 }