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.math.analysis.solvers;
019
020 import org.apache.commons.math.util.FastMath;
021
022 /**
023 * Implements the <em>Secant</em> method for root-finding (approximating a
024 * zero of a univariate real function). The solution that is maintained is
025 * not bracketed, and as such convergence is not guaranteed.
026 *
027 * <p>Implementation based on the following article: M. Dowell and P. Jarratt,
028 * <em>A modified regula falsi method for computing the root of an
029 * equation</em>, BIT Numerical Mathematics, volume 11, number 2,
030 * pages 168-174, Springer, 1971.</p>
031 *
032 * <p>Note that since release 3.0 this class implements the actual
033 * <em>Secant</em> algorithm, and not a modified one. As such, the 3.0 version
034 * is not backwards compatible with previous versions. To use an algorithm
035 * similar to the pre-3.0 releases, use the
036 * {@link IllinoisSolver <em>Illinois</em>} algorithm or the
037 * {@link PegasusSolver <em>Pegasus</em>} algorithm.</p>
038 *
039 * @version $Id: SecantSolver.java 1145146 2011-07-11 12:32:14Z erans $
040 */
041 public class SecantSolver extends AbstractUnivariateRealSolver {
042
043 /** Default absolute accuracy. */
044 protected static final double DEFAULT_ABSOLUTE_ACCURACY = 1e-6;
045
046 /** Construct a solver with default accuracy (1e-6). */
047 public SecantSolver() {
048 super(DEFAULT_ABSOLUTE_ACCURACY);
049 }
050
051 /**
052 * Construct a solver.
053 *
054 * @param absoluteAccuracy absolute accuracy
055 */
056 public SecantSolver(final double absoluteAccuracy) {
057 super(absoluteAccuracy);
058 }
059
060 /**
061 * Construct a solver.
062 *
063 * @param relativeAccuracy relative accuracy
064 * @param absoluteAccuracy absolute accuracy
065 */
066 public SecantSolver(final double relativeAccuracy,
067 final double absoluteAccuracy) {
068 super(relativeAccuracy, absoluteAccuracy);
069 }
070
071 /** {@inheritDoc} */
072 @Override
073 protected final double doSolve() {
074 // Get initial solution
075 double x0 = getMin();
076 double x1 = getMax();
077 double f0 = computeObjectiveValue(x0);
078 double f1 = computeObjectiveValue(x1);
079
080 // If one of the bounds is the exact root, return it. Since these are
081 // not under-approximations or over-approximations, we can return them
082 // regardless of the allowed solutions.
083 if (f0 == 0.0) {
084 return x0;
085 }
086 if (f1 == 0.0) {
087 return x1;
088 }
089
090 // Verify bracketing of initial solution.
091 verifyBracketing(x0, x1);
092
093 // Get accuracies.
094 final double ftol = getFunctionValueAccuracy();
095 final double atol = getAbsoluteAccuracy();
096 final double rtol = getRelativeAccuracy();
097
098 // Keep finding better approximations.
099 while (true) {
100 // Calculate the next approximation.
101 final double x = x1 - ((f1 * (x1 - x0)) / (f1 - f0));
102 final double fx = computeObjectiveValue(x);
103
104 // If the new approximation is the exact root, return it. Since
105 // this is not an under-approximation or an over-approximation,
106 // we can return it regardless of the allowed solutions.
107 if (fx == 0.0) {
108 return x;
109 }
110
111 // Update the bounds with the new approximation.
112 x0 = x1;
113 f0 = f1;
114 x1 = x;
115 f1 = fx;
116
117 // If the function value of the last approximation is too small,
118 // given the function value accuracy, then we can't get closer to
119 // the root than we already are.
120 if (FastMath.abs(f1) <= ftol) {
121 return x1;
122 }
123
124 // If the current interval is within the given accuracies, we
125 // are satisfied with the current approximation.
126 if (FastMath.abs(x1 - x0) < FastMath.max(rtol * FastMath.abs(x1), atol)) {
127 return x1;
128 }
129 }
130 }
131
132 }