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.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 }