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.  *      http://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. package org.apache.commons.math4.legacy.analysis.solvers;

  20. import org.apache.commons.math4.legacy.analysis.UnivariateFunction;
  21. import org.apache.commons.math4.legacy.exception.ConvergenceException;
  22. import org.apache.commons.math4.legacy.exception.MathInternalError;
  23. import org.apache.commons.math4.core.jdkmath.JdkMath;

  25. /**
  26.  * Base class for all bracketing <em>Secant</em>-based methods for root-finding
  27.  * (approximating a zero of a univariate real function).
  28.  *
  29.  * <p>Implementation of the {@link RegulaFalsiSolver <em>Regula Falsi</em>} and
  30.  * {@link IllinoisSolver <em>Illinois</em>} methods is based on the
  31.  * following article: M. Dowell and P. Jarratt,
  32.  * <em>A modified regula falsi method for computing the root of an
  33.  * equation</em>, BIT Numerical Mathematics, volume 11, number 2,
  34.  * pages 168-174, Springer, 1971.</p>
  35.  *
  36.  * <p>Implementation of the {@link PegasusSolver <em>Pegasus</em>} method is
  37.  * based on the following article: M. Dowell and P. Jarratt,
  38.  * <em>The "Pegasus" method for computing the root of an equation</em>,
  39.  * BIT Numerical Mathematics, volume 12, number 4, pages 503-508, Springer,
  40.  * 1972.</p>
  41.  *
  42.  * <p>The {@link SecantSolver <em>Secant</em>} method is <em>not</em> a
  43.  * bracketing method, so it is not implemented here. It has a separate
  44.  * implementation.</p>
  45.  *
  46.  * @since 3.0
  47.  */
  48. public abstract class BaseSecantSolver
  49.     extends AbstractUnivariateSolver
  50.     implements BracketedUnivariateSolver<UnivariateFunction> {

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

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

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

  61.     /**
  62.      * Construct a solver.
  63.      *
  64.      * @param absoluteAccuracy Absolute accuracy.
  65.      * @param method <em>Secant</em>-based root-finding method to use.
  66.      */
  67.     protected BaseSecantSolver(final double absoluteAccuracy, final Method method) {
  68.         super(absoluteAccuracy);
  69.         this.allowed = AllowedSolution.ANY_SIDE;
  70.         this.method = method;
  71.     }

  73.     /**
  74.      * Construct a solver.
  75.      *
  76.      * @param relativeAccuracy Relative accuracy.
  77.      * @param absoluteAccuracy Absolute accuracy.
  78.      * @param method <em>Secant</em>-based root-finding method to use.
  79.      */
  80.     protected BaseSecantSolver(final double relativeAccuracy,
  81.                                final double absoluteAccuracy,
  82.                                final Method method) {
  83.         super(relativeAccuracy, absoluteAccuracy);
  84.         this.allowed = AllowedSolution.ANY_SIDE;
  85.         this.method = method;
  86.     }

  88.     /**
  89.      * Construct a solver.
  90.      *
  91.      * @param relativeAccuracy Maximum relative error.
  92.      * @param absoluteAccuracy Maximum absolute error.
  93.      * @param functionValueAccuracy Maximum function value error.
  94.      * @param method <em>Secant</em>-based root-finding method to use
  95.      */
  96.     protected BaseSecantSolver(final double relativeAccuracy,
  97.                                final double absoluteAccuracy,
  98.                                final double functionValueAccuracy,
  99.                                final Method method) {
  100.         super(relativeAccuracy, absoluteAccuracy, functionValueAccuracy);
  101.         this.allowed = AllowedSolution.ANY_SIDE;
  102.         this.method = method;
  103.     }

  105.     /** {@inheritDoc} */
  106.     @Override
  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.     }

  113.     /** {@inheritDoc} */
  114.     @Override
  115.     public double solve(final int maxEval, final UnivariateFunction f,
  116.                         final double min, final double max, final double startValue,
  117.                         final AllowedSolution allowedSolution) {
  118.         this.allowed = allowedSolution;
  119.         return super.solve(maxEval, f, min, max, startValue);
  120.     }

  122.     /** {@inheritDoc} */
  123.     @Override
  124.     public double solve(final int maxEval, final UnivariateFunction f,
  125.                         final double min, final double max, final double startValue) {
  126.         return solve(maxEval, f, min, max, startValue, AllowedSolution.ANY_SIDE);
  127.     }

  129.     /**
  130.      * {@inheritDoc}
  131.      *
  132.      * @throws ConvergenceException if the algorithm failed due to finite
  133.      * precision.
  134.      */
  135.     @Override
  136.     protected final double doSolve()
  137.         throws ConvergenceException {
  138.         // Get initial solution
  139.         double x0 = getMin();
  140.         double x1 = getMax();
  141.         double f0 = computeObjectiveValue(x0);
  142.         double f1 = computeObjectiveValue(x1);

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

  154.         // Verify bracketing of initial solution.
  155.         verifyBracketing(x0, x1);

  157.         // Get accuracies.
  158.         final double ftol = getFunctionValueAccuracy();
  159.         final double atol = getAbsoluteAccuracy();
  160.         final double rtol = getRelativeAccuracy();

  162.         // Keep track of inverted intervals, meaning that the left bound is
  163.         // larger than the right bound.
  164.         boolean inverted = false;

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

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

  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;

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

  242.             // If the current interval is within the given accuracies, we
  243.             // are satisfied with the current approximation.
  244.             if (JdkMath.abs(x1 - x0) < JdkMath.max(rtol * JdkMath.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.     }

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

  267.         /**
  268.          * The {@link RegulaFalsiSolver <em>Regula Falsi</em>} or
  269.          * <em>False Position</em> method.
  270.          */
  271.         REGULA_FALSI,

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

  276.         /** The {@link PegasusSolver <em>Pegasus</em>} method. */
  277.         PEGASUS;
  278.     }
  279. }
Created with JaCoCo 0.8.10.202304240956Code Coverage Report for Apache Commons Math 4.0-SNAPSHOT