View Javadoc
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 }