BaseSecantSolver.java

  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.  */

  18. /*
  19.  * This is not the original file distributed by the Apache Software Foundation
  20.  * It has been modified by the Hipparchus project
  21.  */

  23. package org.hipparchus.analysis.solvers;

  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;

  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> {

  58.     /** Default absolute accuracy. */
  59.     protected static final double DEFAULT_ABSOLUTE_ACCURACY = 1e-6;

  61.     /** The kinds of solutions that the algorithm may accept. */
  62.     private AllowedSolution allowed;

  64.     /** The <em>Secant</em>-based root-finding method to use. */
  65.     private final Method method;

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

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

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

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

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

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

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

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

  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);

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

  182.         // Verify bracketing of initial solution.
  183.         verifyBracketing(x0, x1);

  185.         // Get accuracies.
  186.         final double ftol = getFunctionValueAccuracy();
  187.         final double atol = getAbsoluteAccuracy();
  188.         final double rtol = getRelativeAccuracy();

  190.         // Keep track of inverted intervals, meaning that the left bound is
  191.         // larger than the right bound.
  192.         boolean inverted = false;

  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);

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

  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;

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

  255.     /** <em>Secant</em>-based root-finding methods. */
  256.     protected enum Method {

  258.         /**
  259.          * The {@link RegulaFalsiSolver <em>Regula Falsi</em>} or
  260.          * <em>False Position</em> method.
  261.          */
  262.         REGULA_FALSI,

  264.         /** The {@link IllinoisSolver <em>Illinois</em>} method. */
  265.         ILLINOIS,

  267.         /** The {@link PegasusSolver <em>Pegasus</em>} method. */
  268.         PEGASUS

  270.     }
  271. }
Created with JaCoCo 0.8.12.202403310830