UnivariateSolverUtils.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.  */
  22. package org.hipparchus.analysis.solvers;

  24. import org.hipparchus.CalculusFieldElement;
  25. import org.hipparchus.analysis.CalculusFieldUnivariateFunction;
  26. import org.hipparchus.analysis.UnivariateFunction;
  27. import org.hipparchus.exception.LocalizedCoreFormats;
  28. import org.hipparchus.exception.MathIllegalArgumentException;
  29. import org.hipparchus.exception.NullArgumentException;
  30. import org.hipparchus.util.FastMath;
  31. import org.hipparchus.util.MathArrays;
  32. import org.hipparchus.util.MathUtils;

  34. /**
  35.  * Utility routines for {@link UnivariateSolver} objects.
  36.  *
  37.  */
  38. public class UnivariateSolverUtils {
  39.     /**
  40.      * Class contains only static methods.
  41.      */
  42.     private UnivariateSolverUtils() {}

  44.     /**
  45.      * Convenience method to find a zero of a univariate real function.  A default
  46.      * solver is used.
  47.      *
  48.      * @param function Function.
  49.      * @param x0 Lower bound for the interval.
  50.      * @param x1 Upper bound for the interval.
  51.      * @return a value where the function is zero.
  52.      * @throws MathIllegalArgumentException if the function has the same sign at the
  53.      * endpoints.
  54.      * @throws NullArgumentException if {@code function} is {@code null}.
  55.      */
  56.     public static double solve(UnivariateFunction function, double x0, double x1)
  57.         throws MathIllegalArgumentException, NullArgumentException {
  58.         MathUtils.checkNotNull(function, LocalizedCoreFormats.FUNCTION);
  59.         final UnivariateSolver solver = new BrentSolver();
  60.         return solver.solve(Integer.MAX_VALUE, function, x0, x1);
  61.     }

  63.     /**
  64.      * Convenience method to find a zero of a univariate real function.  A default
  65.      * solver is used.
  66.      *
  67.      * @param function Function.
  68.      * @param x0 Lower bound for the interval.
  69.      * @param x1 Upper bound for the interval.
  70.      * @param absoluteAccuracy Accuracy to be used by the solver.
  71.      * @return a value where the function is zero.
  72.      * @throws MathIllegalArgumentException if the function has the same sign at the
  73.      * endpoints.
  74.      * @throws NullArgumentException if {@code function} is {@code null}.
  75.      */
  76.     public static double solve(UnivariateFunction function,
  77.                                double x0, double x1,
  78.                                double absoluteAccuracy)
  79.         throws MathIllegalArgumentException, NullArgumentException {
  80.         MathUtils.checkNotNull(function, LocalizedCoreFormats.FUNCTION);
  81.         final UnivariateSolver solver = new BrentSolver(absoluteAccuracy);
  82.         return solver.solve(Integer.MAX_VALUE, function, x0, x1);
  83.     }

  85.     /**
  86.      * Force a root found by a non-bracketing solver to lie on a specified side,
  87.      * as if the solver were a bracketing one.
  88.      *
  89.      * @param maxEval maximal number of new evaluations of the function
  90.      * (evaluations already done for finding the root should have already been subtracted
  91.      * from this number)
  92.      * @param f function to solve
  93.      * @param bracketing bracketing solver to use for shifting the root
  94.      * @param baseRoot original root found by a previous non-bracketing solver
  95.      * @param min minimal bound of the search interval
  96.      * @param max maximal bound of the search interval
  97.      * @param allowedSolution the kind of solutions that the root-finding algorithm may
  98.      * accept as solutions.
  99.      * @return a root approximation, on the specified side of the exact root
  100.      * @throws MathIllegalArgumentException if the function has the same sign at the
  101.      * endpoints.
  102.      */
  103.     public static double forceSide(final int maxEval, final UnivariateFunction f,
  104.                                    final BracketedUnivariateSolver<UnivariateFunction> bracketing,
  105.                                    final double baseRoot, final double min, final double max,
  106.                                    final AllowedSolution allowedSolution)
  107.         throws MathIllegalArgumentException {

  109.         if (allowedSolution == AllowedSolution.ANY_SIDE) {
  110.             // no further bracketing required
  111.             return baseRoot;
  112.         }

  114.         // find a very small interval bracketing the root
  115.         final double step = FastMath.max(bracketing.getAbsoluteAccuracy(),
  116.                                          FastMath.abs(baseRoot * bracketing.getRelativeAccuracy()));
  117.         double xLo        = FastMath.max(min, baseRoot - step);
  118.         double fLo        = f.value(xLo);
  119.         double xHi        = FastMath.min(max, baseRoot + step);
  120.         double fHi        = f.value(xHi);
  121.         int remainingEval = maxEval - 2;
  122.         while (remainingEval > 0) {

  124.             if ((fLo >= 0 && fHi <= 0) || (fLo <= 0 && fHi >= 0)) {
  125.                 // compute the root on the selected side
  126.                 return bracketing.solve(remainingEval, f, xLo, xHi, baseRoot, allowedSolution);
  127.             }

  129.             // try increasing the interval
  130.             boolean changeLo = false;
  131.             boolean changeHi = false;
  132.             if (fLo < fHi) {
  133.                 // increasing function
  134.                 if (fLo >= 0) {
  135.                     changeLo = true;
  136.                 } else {
  137.                     changeHi = true;
  138.                 }
  139.             } else if (fLo > fHi) {
  140.                 // decreasing function
  141.                 if (fLo <= 0) {
  142.                     changeLo = true;
  143.                 } else {
  144.                     changeHi = true;
  145.                 }
  146.             } else {
  147.                 // unknown variation
  148.                 changeLo = true;
  149.                 changeHi = true;
  150.             }

  152.             // update the lower bound
  153.             if (changeLo) {
  154.                 xLo = FastMath.max(min, xLo - step);
  155.                 fLo  = f.value(xLo);
  156.                 remainingEval--;
  157.             }

  159.             // update the higher bound
  160.             if (changeHi) {
  161.                 xHi = FastMath.min(max, xHi + step);
  162.                 fHi  = f.value(xHi);
  163.                 remainingEval--;
  164.             }

  166.         }

  168.         throw new MathIllegalArgumentException(LocalizedCoreFormats.FAILED_BRACKETING,
  169.                                                xLo, xHi, fLo, fHi,
  170.                                                maxEval - remainingEval, maxEval, baseRoot,
  171.                                                min, max);

  173.     }

  175.     /**
  176.      * This method simply calls {@link #bracket(UnivariateFunction, double, double, double,
  177.      * double, double, int) bracket(function, initial, lowerBound, upperBound, q, r, maximumIterations)}
  178.      * with {@code q} and {@code r} set to 1.0 and {@code maximumIterations} set to {@code Integer.MAX_VALUE}.
  179.      * <p>
  180.      * <strong>Note: </strong> this method can take {@code Integer.MAX_VALUE}
  181.      * iterations to throw a {@code MathIllegalStateException.}  Unless you are
  182.      * confident that there is a root between {@code lowerBound} and
  183.      * {@code upperBound} near {@code initial}, it is better to use
  184.      * {@link #bracket(UnivariateFunction, double, double, double, double,double, int)
  185.      * bracket(function, initial, lowerBound, upperBound, q, r, maximumIterations)},
  186.      * explicitly specifying the maximum number of iterations.</p>
  187.      *
  188.      * @param function Function.
  189.      * @param initial Initial midpoint of interval being expanded to
  190.      * bracket a root.
  191.      * @param lowerBound Lower bound (a is never lower than this value)
  192.      * @param upperBound Upper bound (b never is greater than this
  193.      * value).
  194.      * @return a two-element array holding a and b.
  195.      * @throws MathIllegalArgumentException if a root cannot be bracketed.
  196.      * @throws MathIllegalArgumentException if {@code maximumIterations <= 0}.
  197.      * @throws NullArgumentException if {@code function} is {@code null}.
  198.      */
  199.     public static double[] bracket(UnivariateFunction function,
  200.                                    double initial,
  201.                                    double lowerBound, double upperBound)
  202.         throws MathIllegalArgumentException, NullArgumentException {
  203.         return bracket(function, initial, lowerBound, upperBound, 1.0, 1.0, Integer.MAX_VALUE);
  204.     }

  206.      /**
  207.      * This method simply calls {@link #bracket(UnivariateFunction, double, double, double,
  208.      * double, double, int) bracket(function, initial, lowerBound, upperBound, q, r, maximumIterations)}
  209.      * with {@code q} and {@code r} set to 1.0.
  210.      * @param function Function.
  211.      * @param initial Initial midpoint of interval being expanded to
  212.      * bracket a root.
  213.      * @param lowerBound Lower bound (a is never lower than this value).
  214.      * @param upperBound Upper bound (b never is greater than this
  215.      * value).
  216.      * @param maximumIterations Maximum number of iterations to perform
  217.      * @return a two element array holding a and b.
  218.      * @throws MathIllegalArgumentException if the algorithm fails to find a and b
  219.      * satisfying the desired conditions.
  220.      * @throws MathIllegalArgumentException if {@code maximumIterations <= 0}.
  221.      * @throws NullArgumentException if {@code function} is {@code null}.
  222.      */
  223.     public static double[] bracket(UnivariateFunction function,
  224.                                    double initial,
  225.                                    double lowerBound, double upperBound,
  226.                                    int maximumIterations)
  227.         throws MathIllegalArgumentException, NullArgumentException {
  228.         return bracket(function, initial, lowerBound, upperBound, 1.0, 1.0, maximumIterations);
  229.     }

  231.     /**
  232.      * This method attempts to find two values a and b satisfying <ul>
  233.      * <li> {@code lowerBound <= a < initial < b <= upperBound} </li>
  234.      * <li> {@code f(a) * f(b) <= 0} </li>
  235.      * </ul>
  236.      * If {@code f} is continuous on {@code [a,b]}, this means that {@code a}
  237.      * and {@code b} bracket a root of {@code f}.
  238.      * <p>
  239.      * The algorithm checks the sign of \( f(l_k) \) and \( f(u_k) \) for increasing
  240.      * values of k, where \( l_k = max(lower, initial - \delta_k) \),
  241.      * \( u_k = min(upper, initial + \delta_k) \), using recurrence
  242.      * \( \delta_{k+1} = r \delta_k + q, \delta_0 = 0\) and starting search with \( k=1 \).
  243.      * The algorithm stops when one of the following happens: <ul>
  244.      * <li> at least one positive and one negative value have been found --  success!</li>
  245.      * <li> both endpoints have reached their respective limits -- MathIllegalArgumentException </li>
  246.      * <li> {@code maximumIterations} iterations elapse -- MathIllegalArgumentException </li></ul>
  247.      * <p>
  248.      * If different signs are found at first iteration ({@code k=1}), then the returned
  249.      * interval will be \( [a, b] = [l_1, u_1] \). If different signs are found at a later
  250.      * iteration {@code k>1}, then the returned interval will be either
  251.      * \( [a, b] = [l_{k+1}, l_{k}] \) or \( [a, b] = [u_{k}, u_{k+1}] \). A root solver called
  252.      * with these parameters will therefore start with the smallest bracketing interval known
  253.      * at this step.
  254.      * </p>
  255.      * <p>
  256.      * Interval expansion rate is tuned by changing the recurrence parameters {@code r} and
  257.      * {@code q}. When the multiplicative factor {@code r} is set to 1, the sequence is a
  258.      * simple arithmetic sequence with linear increase. When the multiplicative factor {@code r}
  259.      * is larger than 1, the sequence has an asymptotically exponential rate. Note than the
  260.      * additive parameter {@code q} should never be set to zero, otherwise the interval would
  261.      * degenerate to the single initial point for all values of {@code k}.
  262.      * </p>
  263.      * <p>
  264.      * As a rule of thumb, when the location of the root is expected to be approximately known
  265.      * within some error margin, {@code r} should be set to 1 and {@code q} should be set to the
  266.      * order of magnitude of the error margin. When the location of the root is really a wild guess,
  267.      * then {@code r} should be set to a value larger than 1 (typically 2 to double the interval
  268.      * length at each iteration) and {@code q} should be set according to half the initial
  269.      * search interval length.
  270.      * </p>
  271.      * <p>
  272.      * As an example, if we consider the trivial function {@code f(x) = 1 - x} and use
  273.      * {@code initial = 4}, {@code r = 1}, {@code q = 2}, the algorithm will compute
  274.      * {@code f(4-2) = f(2) = -1} and {@code f(4+2) = f(6) = -5} for {@code k = 1}, then
  275.      * {@code f(4-4) = f(0) = +1} and {@code f(4+4) = f(8) = -7} for {@code k = 2}. Then it will
  276.      * return the interval {@code [0, 2]} as the smallest one known to be bracketing the root.
  277.      * As shown by this example, the initial value (here {@code 4}) may lie outside of the returned
  278.      * bracketing interval.
  279.      * </p>
  280.      * @param function function to check
  281.      * @param initial Initial midpoint of interval being expanded to
  282.      * bracket a root.
  283.      * @param lowerBound Lower bound (a is never lower than this value).
  284.      * @param upperBound Upper bound (b never is greater than this
  285.      * value).
  286.      * @param q additive offset used to compute bounds sequence (must be strictly positive)
  287.      * @param r multiplicative factor used to compute bounds sequence
  288.      * @param maximumIterations Maximum number of iterations to perform
  289.      * @return a two element array holding the bracketing values.
  290.      * @exception MathIllegalArgumentException if function cannot be bracketed in the search interval
  291.      */
  292.     public static double[] bracket(final UnivariateFunction function, final double initial,
  293.                                    final double lowerBound, final double upperBound,
  294.                                    final double q, final double r, final int maximumIterations)
  295.         throws MathIllegalArgumentException {

  297.         MathUtils.checkNotNull(function, LocalizedCoreFormats.FUNCTION);

  299.         if (q <= 0)  {
  300.             throw new MathIllegalArgumentException(LocalizedCoreFormats.NUMBER_TOO_SMALL_BOUND_EXCLUDED,
  301.                                                    q, 0);
  302.         }
  303.         if (maximumIterations <= 0)  {
  304.             throw new MathIllegalArgumentException(LocalizedCoreFormats.INVALID_MAX_ITERATIONS, maximumIterations);
  305.         }
  306.         verifySequence(lowerBound, initial, upperBound);

  308.         // initialize the recurrence
  309.         double a     = initial;
  310.         double b     = initial;
  311.         double fa    = Double.NaN;
  312.         double fb    = Double.NaN;
  313.         double delta = 0;

  315.         for (int numIterations = 0;
  316.              (numIterations < maximumIterations) && (a > lowerBound || b < upperBound);
  317.              ++numIterations) {

  319.             final double previousA  = a;
  320.             final double previousFa = fa;
  321.             final double previousB  = b;
  322.             final double previousFb = fb;

  324.             delta = r * delta + q;
  325.             a     = FastMath.max(initial - delta, lowerBound);
  326.             b     = FastMath.min(initial + delta, upperBound);
  327.             fa    = function.value(a);
  328.             fb    = function.value(b);

  330.             if (numIterations == 0) {
  331.                 // at first iteration, we don't have a previous interval
  332.                 // we simply compare both sides of the initial interval
  333.                 if (fa * fb <= 0) {
  334.                     // the first interval already brackets a root
  335.                     return new double[] { a, b };
  336.                 }
  337.             } else {
  338.                 // we have a previous interval with constant sign and expand it,
  339.                 // we expect sign changes to occur at boundaries
  340.                 if (fa * previousFa <= 0) {
  341.                     // sign change detected at near lower bound
  342.                     return new double[] { a, previousA };
  343.                 } else if (fb * previousFb <= 0) {
  344.                     // sign change detected at near upper bound
  345.                     return new double[] { previousB, b };
  346.                 }
  347.             }

  349.         }

  351.         // no bracketing found
  352.         throw new MathIllegalArgumentException(LocalizedCoreFormats.NOT_BRACKETING_INTERVAL,
  353.                                                a, b, fa, fb);

  355.     }

  357.     /**
  358.      * This method simply calls {@link #bracket(CalculusFieldUnivariateFunction,
  359.      * CalculusFieldElement, CalculusFieldElement, CalculusFieldElement, CalculusFieldElement,
  360.      * CalculusFieldElement, int) bracket(function, initial, lowerBound, upperBound, q, r, maximumIterations)}
  361.      * with {@code q} and {@code r} set to 1.0 and {@code maximumIterations} set to {@code Integer.MAX_VALUE}.
  362.      * <p>
  363.      * <strong>Note: </strong> this method can take {@code Integer.MAX_VALUE}
  364.      * iterations to throw a {@code MathIllegalStateException.}  Unless you are
  365.      * confident that there is a root between {@code lowerBound} and
  366.      * {@code upperBound} near {@code initial}, it is better to use
  367.      * {@link #bracket(UnivariateFunction, double, double, double, double,double, int)
  368.      * bracket(function, initial, lowerBound, upperBound, q, r, maximumIterations)},
  369.      * explicitly specifying the maximum number of iterations.</p>
  370.      *
  371.      * @param function Function.
  372.      * @param initial Initial midpoint of interval being expanded to
  373.      * bracket a root.
  374.      * @param lowerBound Lower bound (a is never lower than this value)
  375.      * @param upperBound Upper bound (b never is greater than this
  376.      * value).
  377.      * @param <T> type of the field elements
  378.      * @return a two-element array holding a and b.
  379.      * @throws MathIllegalArgumentException if a root cannot be bracketed.
  380.      * @throws MathIllegalArgumentException if {@code maximumIterations <= 0}.
  381.      * @throws NullArgumentException if {@code function} is {@code null}.
  382.      * @since 1.2
  383.      */
  384.     public static <T extends CalculusFieldElement<T>> T[] bracket(CalculusFieldUnivariateFunction<T> function,
  385.                                                               T initial,
  386.                                                               T lowerBound, T upperBound)
  387.         throws MathIllegalArgumentException, NullArgumentException {
  388.         return bracket(function, initial, lowerBound, upperBound,
  389.                        initial.getField().getOne(), initial.getField().getOne(),
  390.                        Integer.MAX_VALUE);
  391.     }

  393.      /**
  394.      * This method simply calls {@link #bracket(CalculusFieldUnivariateFunction,
  395.      * CalculusFieldElement, CalculusFieldElement, CalculusFieldElement, CalculusFieldElement,
  396.      * CalculusFieldElement, int) bracket(function, initial, lowerBound, upperBound, q, r, maximumIterations)}
  397.      * with {@code q} and {@code r} set to 1.0.
  398.      * @param function Function.
  399.      * @param initial Initial midpoint of interval being expanded to
  400.      * bracket a root.
  401.      * @param lowerBound Lower bound (a is never lower than this value).
  402.      * @param upperBound Upper bound (b never is greater than this
  403.      * value).
  404.      * @param maximumIterations Maximum number of iterations to perform
  405.      * @param <T> type of the field elements
  406.      * @return a two element array holding a and b.
  407.      * @throws MathIllegalArgumentException if the algorithm fails to find a and b
  408.      * satisfying the desired conditions.
  409.      * @throws MathIllegalArgumentException if {@code maximumIterations <= 0}.
  410.      * @throws NullArgumentException if {@code function} is {@code null}.
  411.      * @since 1.2
  412.      */
  413.     public static <T extends CalculusFieldElement<T>> T[] bracket(CalculusFieldUnivariateFunction<T> function,
  414.                                                               T initial,
  415.                                                               T lowerBound, T upperBound,
  416.                                                               int maximumIterations)
  417.         throws MathIllegalArgumentException, NullArgumentException {
  418.         return bracket(function, initial, lowerBound, upperBound,
  419.                        initial.getField().getOne(), initial.getField().getOne(),
  420.                        maximumIterations);
  421.     }

  423.     /**
  424.      * This method attempts to find two values a and b satisfying <ul>
  425.      * <li> {@code lowerBound <= a < initial < b <= upperBound} </li>
  426.      * <li> {@code f(a) * f(b) <= 0} </li>
  427.      * </ul>
  428.      * If {@code f} is continuous on {@code [a,b]}, this means that {@code a}
  429.      * and {@code b} bracket a root of {@code f}.
  430.      * <p>
  431.      * The algorithm checks the sign of \( f(l_k) \) and \( f(u_k) \) for increasing
  432.      * values of k, where \( l_k = max(lower, initial - \delta_k) \),
  433.      * \( u_k = min(upper, initial + \delta_k) \), using recurrence
  434.      * \( \delta_{k+1} = r \delta_k + q, \delta_0 = 0\) and starting search with \( k=1 \).
  435.      * The algorithm stops when one of the following happens: <ul>
  436.      * <li> at least one positive and one negative value have been found --  success!</li>
  437.      * <li> both endpoints have reached their respective limits -- MathIllegalArgumentException </li>
  438.      * <li> {@code maximumIterations} iterations elapse -- MathIllegalArgumentException </li></ul>
  439.      * <p>
  440.      * If different signs are found at first iteration ({@code k=1}), then the returned
  441.      * interval will be \( [a, b] = [l_1, u_1] \). If different signs are found at a later
  442.      * iteration {@code k>1}, then the returned interval will be either
  443.      * \( [a, b] = [l_{k+1}, l_{k}] \) or \( [a, b] = [u_{k}, u_{k+1}] \). A root solver called
  444.      * with these parameters will therefore start with the smallest bracketing interval known
  445.      * at this step.
  446.      * </p>
  447.      * <p>
  448.      * Interval expansion rate is tuned by changing the recurrence parameters {@code r} and
  449.      * {@code q}. When the multiplicative factor {@code r} is set to 1, the sequence is a
  450.      * simple arithmetic sequence with linear increase. When the multiplicative factor {@code r}
  451.      * is larger than 1, the sequence has an asymptotically exponential rate. Note than the
  452.      * additive parameter {@code q} should never be set to zero, otherwise the interval would
  453.      * degenerate to the single initial point for all values of {@code k}.
  454.      * </p>
  455.      * <p>
  456.      * As a rule of thumb, when the location of the root is expected to be approximately known
  457.      * within some error margin, {@code r} should be set to 1 and {@code q} should be set to the
  458.      * order of magnitude of the error margin. When the location of the root is really a wild guess,
  459.      * then {@code r} should be set to a value larger than 1 (typically 2 to double the interval
  460.      * length at each iteration) and {@code q} should be set according to half the initial
  461.      * search interval length.
  462.      * </p>
  463.      * <p>
  464.      * As an example, if we consider the trivial function {@code f(x) = 1 - x} and use
  465.      * {@code initial = 4}, {@code r = 1}, {@code q = 2}, the algorithm will compute
  466.      * {@code f(4-2) = f(2) = -1} and {@code f(4+2) = f(6) = -5} for {@code k = 1}, then
  467.      * {@code f(4-4) = f(0) = +1} and {@code f(4+4) = f(8) = -7} for {@code k = 2}. Then it will
  468.      * return the interval {@code [0, 2]} as the smallest one known to be bracketing the root.
  469.      * As shown by this example, the initial value (here {@code 4}) may lie outside of the returned
  470.      * bracketing interval.
  471.      * </p>
  472.      * @param function function to check
  473.      * @param initial Initial midpoint of interval being expanded to
  474.      * bracket a root.
  475.      * @param lowerBound Lower bound (a is never lower than this value).
  476.      * @param upperBound Upper bound (b never is greater than this
  477.      * value).
  478.      * @param q additive offset used to compute bounds sequence (must be strictly positive)
  479.      * @param r multiplicative factor used to compute bounds sequence
  480.      * @param maximumIterations Maximum number of iterations to perform
  481.      * @param <T> type of the field elements
  482.      * @return a two element array holding the bracketing values.
  483.      * @exception MathIllegalArgumentException if function cannot be bracketed in the search interval
  484.      * @since 1.2
  485.      */
  486.     public static <T extends CalculusFieldElement<T>> T[] bracket(final CalculusFieldUnivariateFunction<T> function,
  487.                                                               final T initial,
  488.                                                               final T lowerBound, final T upperBound,
  489.                                                               final T q, final T r,
  490.                                                               final int maximumIterations)
  491.         throws MathIllegalArgumentException {

  493.         MathUtils.checkNotNull(function, LocalizedCoreFormats.FUNCTION);

  495.         if (q.getReal() <= 0)  {
  496.             throw new MathIllegalArgumentException(LocalizedCoreFormats.NUMBER_TOO_SMALL_BOUND_EXCLUDED,
  497.                                                    q, 0);
  498.         }
  499.         if (maximumIterations <= 0)  {
  500.             throw new MathIllegalArgumentException(LocalizedCoreFormats.INVALID_MAX_ITERATIONS, maximumIterations);
  501.         }
  502.         verifySequence(lowerBound.getReal(), initial.getReal(), upperBound.getReal());

  504.         // initialize the recurrence
  505.         T a     = initial;
  506.         T b     = initial;
  507.         T fa    = null;
  508.         T fb    = null;
  509.         T delta = initial.getField().getZero();

  511.         for (int numIterations = 0;
  512.              (numIterations < maximumIterations) &&
  513.              (a.getReal() > lowerBound.getReal() || b.getReal() < upperBound.getReal());
  514.              ++numIterations) {

  516.             final T previousA  = a;
  517.             final T previousFa = fa;
  518.             final T previousB  = b;
  519.             final T previousFb = fb;

  521.             delta = r.multiply(delta).add(q);
  522.             a     = max(initial.subtract(delta), lowerBound);
  523.             b     = min(initial.add(delta), upperBound);
  524.             fa    = function.value(a);
  525.             fb    = function.value(b);

  527.             if (numIterations == 0) {
  528.                 // at first iteration, we don't have a previous interval
  529.                 // we simply compare both sides of the initial interval
  530.                 if (fa.multiply(fb).getReal() <= 0) {
  531.                     // the first interval already brackets a root
  532.                     final T[] interval = MathArrays.buildArray(initial.getField(), 2);
  533.                     interval[0] = a;
  534.                     interval[1] = b;
  535.                     return interval;
  536.                 }
  537.             } else {
  538.                 // we have a previous interval with constant sign and expand it,
  539.                 // we expect sign changes to occur at boundaries
  540.                 if (fa.multiply(previousFa).getReal() <= 0) {
  541.                     // sign change detected at near lower bound
  542.                     final T[] interval = MathArrays.buildArray(initial.getField(), 2);
  543.                     interval[0] = a;
  544.                     interval[1] = previousA;
  545.                     return interval;
  546.                 } else if (fb.multiply(previousFb).getReal() <= 0) {
  547.                     // sign change detected at near upper bound
  548.                     final T[] interval = MathArrays.buildArray(initial.getField(), 2);
  549.                     interval[0] = previousB;
  550.                     interval[1] = b;
  551.                     return interval;
  552.                 }
  553.             }

  555.         }

  557.         // no bracketing found
  558.         throw new MathIllegalArgumentException(LocalizedCoreFormats.NOT_BRACKETING_INTERVAL,
  559.                                                a.getReal(), b.getReal(), fa.getReal(), fb.getReal());

  561.     }

  563.     /** Compute the maximum of two values
  564.      * @param a first value
  565.      * @param b second value
  566.      * @param <T> type of the field elements
  567.      * @return b if a is lesser or equal to b, a otherwise
  568.      * @since 1.2
  569.      */
  570.     private static <T extends CalculusFieldElement<T>> T max(final T a, final T b) {
  571.         return (a.subtract(b).getReal() <= 0) ? b : a;
  572.     }

  574.     /** Compute the minimum of two values
  575.      * @param a first value
  576.      * @param b second value
  577.      * @param <T> type of the field elements
  578.      * @return a if a is lesser or equal to b, b otherwise
  579.      * @since 1.2
  580.      */
  581.     private static <T extends CalculusFieldElement<T>> T min(final T a, final T b) {
  582.         return (a.subtract(b).getReal() <= 0) ? a : b;
  583.     }

  585.     /**
  586.      * Compute the midpoint of two values.
  587.      *
  588.      * @param a first value.
  589.      * @param b second value.
  590.      * @return the midpoint.
  591.      */
  592.     public static double midpoint(double a, double b) {
  593.         return (a + b) * 0.5;
  594.     }

  596.     /**
  597.      * Check whether the interval bounds bracket a root. That is, if the
  598.      * values at the endpoints are not equal to zero, then the function takes
  599.      * opposite signs at the endpoints.
  600.      *
  601.      * @param function Function.
  602.      * @param lower Lower endpoint.
  603.      * @param upper Upper endpoint.
  604.      * @return {@code true} if the function values have opposite signs at the
  605.      * given points.
  606.      * @throws NullArgumentException if {@code function} is {@code null}.
  607.      */
  608.     public static boolean isBracketing(UnivariateFunction function,
  609.                                        final double lower,
  610.                                        final double upper)
  611.         throws NullArgumentException {
  612.         MathUtils.checkNotNull(function, LocalizedCoreFormats.FUNCTION);
  613.         final double fLo = function.value(lower);
  614.         final double fHi = function.value(upper);
  615.         return (fLo >= 0 && fHi <= 0) || (fLo <= 0 && fHi >= 0);
  616.     }

  618.     /**
  619.      * Check whether the arguments form a (strictly) increasing sequence.
  620.      *
  621.      * @param start First number.
  622.      * @param mid Second number.
  623.      * @param end Third number.
  624.      * @return {@code true} if the arguments form an increasing sequence.
  625.      */
  626.     public static boolean isSequence(final double start,
  627.                                      final double mid,
  628.                                      final double end) {
  629.         return (start < mid) && (mid < end);
  630.     }

  632.     /**
  633.      * Check that the endpoints specify an interval.
  634.      *
  635.      * @param lower Lower endpoint.
  636.      * @param upper Upper endpoint.
  637.      * @throws MathIllegalArgumentException if {@code lower >= upper}.
  638.      */
  639.     public static void verifyInterval(final double lower,
  640.                                       final double upper)
  641.         throws MathIllegalArgumentException {
  642.         if (lower >= upper) {
  643.             throw new MathIllegalArgumentException(LocalizedCoreFormats.ENDPOINTS_NOT_AN_INTERVAL,
  644.                                                 lower, upper, false);
  645.         }
  646.     }

  648.     /**
  649.      * Check that {@code lower < initial < upper}.
  650.      *
  651.      * @param lower Lower endpoint.
  652.      * @param initial Initial value.
  653.      * @param upper Upper endpoint.
  654.      * @throws MathIllegalArgumentException if {@code lower >= initial} or
  655.      * {@code initial >= upper}.
  656.      */
  657.     public static void verifySequence(final double lower,
  658.                                       final double initial,
  659.                                       final double upper)
  660.         throws MathIllegalArgumentException {
  661.         verifyInterval(lower, initial);
  662.         verifyInterval(initial, upper);
  663.     }

  665.     /**
  666.      * Check that the endpoints specify an interval and the end points
  667.      * bracket a root.
  668.      *
  669.      * @param function Function.
  670.      * @param lower Lower endpoint.
  671.      * @param upper Upper endpoint.
  672.      * @throws MathIllegalArgumentException if the function has the same sign at the
  673.      * endpoints.
  674.      * @throws NullArgumentException if {@code function} is {@code null}.
  675.      */
  676.     public static void verifyBracketing(UnivariateFunction function,
  677.                                         final double lower,
  678.                                         final double upper)
  679.         throws MathIllegalArgumentException, NullArgumentException {
  680.         MathUtils.checkNotNull(function, LocalizedCoreFormats.FUNCTION);
  681.         verifyInterval(lower, upper);
  682.         if (!isBracketing(function, lower, upper)) {
  683.             throw new MathIllegalArgumentException(LocalizedCoreFormats.NOT_BRACKETING_INTERVAL,
  684.                                                    lower, upper,
  685.                                                    function.value(lower), function.value(upper));
  686.         }
  687.     }
  688. }
Created with JaCoCo 0.8.12.202403310830