Java- Commons.Math BrentSolver returns NoBracketingException。这是一个错误吗？

Question

我写了一个测试程序来复习BrentSolver class through the Apache Commons Math library。

import java.util.TreeSet;
import org.apache.commons.math3.analysis.UnivariateFunction;
import org.apache.commons.math3.analysis.solvers.*;

public class TestBrent {
    public static void main(String[] args) {
        BrentSolver test2 = new BrentSolver(1E-10);       
        UnivariateFunction func = (double x) -> Math.sin(x);

        TreeSet<Double> set = new TreeSet<>();
        for (int i = 1; i <= 500; i++) {
            set.add(test2.solve(1000, func, i, i+1));
        }
        for (Double s : set) {
            if(s > 0)
            System.out.println(s);
        }
    }
}

当运行程序时，返回如下错误

Exception in thread "main" org.apache.commons.math3.exception.NoBracketingException: function values at endpoints do not have different signs, endpoints: [1, 2], values: [0.841, 0.909]
    at org.apache.commons.math3.analysis.solvers.BrentSolver.doSolve(BrentSolver.java:133)
    at org.apache.commons.math3.analysis.solvers.BaseAbstractUnivariateSolver.solve(BaseAbstractUnivariateSolver.java:199)
    at org.apache.commons.math3.analysis.solvers.BaseAbstractUnivariateSolver.solve(BaseAbstractUnivariateSolver.java:204)
    at TestBrent.main(TestBrent.java:12)
Java Result: 1

删除 if 语句允许程序找到根。

import java.util.TreeSet;
import org.apache.commons.math3.analysis.UnivariateFunction;
import org.apache.commons.math3.analysis.solvers.*;

public class TestBrent {
    public static void main(String[] args) {
        BrentSolver test2 = new BrentSolver(1E-10);       
        UnivariateFunction func = (double x) -> Math.sin(x);

        System.out.println(test2.solve(1000, func, 1, 4));

    }
}

这是 Apache Commons 数学库中的错误吗？所有求根算法（牛顿法之外的算法）似乎都有相同的错误。

Answer 1

不，这不是错误。该错误告诉您您的函数（在本例中为 sin(x)）在 x==1 处计算的符号与在 x==2 处计算的符号相同。这表明您传递给它的间隔在其中没有根。为了可靠地工作，求解器需要传递一个包含一个根的区间。对此的粗略检查是区间边界函数的符号不同（当然，这并不排除包含多个根的区间......）。见 the docs，它说：

A solver may require that the interval brackets a single zero root.

在您的第一个示例中，for 循环将传递求解器不同的间隔（[1,2]、[2,3] 等等）。有些人会扎根，有些则不会。不幸的是，for 循环永远不会超过第一次迭代，因为 1 和 2 之间没有根。

在您的第二个示例中，您将区间 [1,4] 传递给求解器。当然，sin(x) 在 x=PI.

的那些点之间确实有一个根

导致错误的是你的问题条件，而不是库中的错误。

如果您想使用第一个公式 - 您可以将行包围起来

set.add(test2.solve(1000, func, i, i+1));

使用 try...catch 循环，例如：

try {
    set.add(test2.solve(1000, func, i, i+1));
}
catch (NoBracketingException e) {
    System.err.println("Oops - error.  Likely cause: No root in interval ["+i+","+(i+1)+"]");
    System.err.println(e.message());
}

Answer 2

这是解决此问题的另一种方法...

/**************************************************************************
**
**    Riemann-Siegel Formula for roots of Zeta(s) on critical line.
**
**************************************************************************
**    Axion004
**    07/31/2015
**
**    This program finds the roots of Zeta(s) using the well known Riemann-
**    Siegel formula. The Riemann–Siegel theta function is approximated 
**    using Stirling's approximation. It also uses an interpolation method to
**    locate zeroes. The coefficients for R(t) are handled by the Taylor
**    Series approximation originally listed by Haselgrove in 1960. It is 
**    necessary to use these coefficients in order to increase computational 
**    speed.
**************************************************************************/

import java.util.TreeSet;
//These two imports are from the Apache Commons Math library
import org.apache.commons.math3.analysis.UnivariateFunction;
import org.apache.commons.math3.analysis.solvers.BracketingNthOrderBrentSolver;

public class RiemannSiegel{
    public static void main(String[] args){
        SiegelMain();
    }

    // Main method
    public static void SiegelMain() {
        System.out.println("Zeroes inside the critical line for " +
                "Zeta(1/2 + it). The t values are referenced below.");
        System.out.println();
        findRoots();
    }

    /**
     * The sign of a calculated double value.
     * @param x - the double value.
     * @return the sign in -1,  1, or 0 format.
    */
    private static int sign(double x) {
    if (x < 0.0)
            return -1;
        else if (x > 0.0)
            return 1;
        else
            return 0;
    }

    /**
     * Finds the roots of a specified function through the 
         * BracketingNthOrderBrentSolver from the Apache Commons Math library.
         * See http://commons.apache.org/proper/commons-math/apidocs/org/
         * apache/commons/math3/analysis/solvers/BracketingNthOrderBrentSolver
         * .html
     * The zeroes inside the interval of 0 < t < 10000 are printed from
         * a TreeSet.
    */
    public static void findRoots() {
    BracketingNthOrderBrentSolver f = new 
            BracketingNthOrderBrentSolver();
        UnivariateFunction func = (double x) -> RiemennZ(x, 4);
        TreeSet<Double> set = new TreeSet<>();
        double i = 1.0;
        while (i < 10000) {
            i+= 0.1;
            if(sign(func.value(i)) != sign(func.value(i+0.1))) {
            set.add(f.solve(1000, func, i, i+0.1));
            }
        }
        set.stream().filter((s) -> (s > 0)).forEach((s) -> {
            System.out.println(s);
        });
    }

    /**
     * Riemann-Siegel theta function using the approximation by the 
         * Stirling series.
     * @param t - the value of t inside the Z(t) function.
     * @return Stirling's approximation for theta(t).
    */
    public static double theta (double t) {
        return (t/2.0 * Math.log(t/(2.0*Math.PI)) - t/2.0 - Math.PI/8.0
                + 1.0/(48.0*Math.pow(t, 1)) + 7.0/(5760*Math.pow(t, 3)));
    }

    /**
     * Computes Math.Floor of the absolute value term passed in as t.
     * @param t - the value of t inside the Z(t) function.
     * @return Math.floor of the absolute value of t.
    */
    public static double fAbs(double t) {
        return Math.floor(Math.abs(t));

    }

    /**
     * Riemann-Siegel Z(t) function implemented per the Riemenn Siegel 
         * formula. See http://mathworld.wolfram.com/Riemann-SiegelFormula.html 
         * for details
     * @param t - the value of t inside the Z(t) function.
         * @param r - referenced for calculating the remainder terms by the
         * Taylor series approximations.
     * @return the approximate value of Z(t) through the Riemann-Siegel
         * formula
    */
    public static double RiemennZ(double t, int r) {

        double twopi = Math.PI * 2.0; 
        double val = Math.sqrt(t/twopi);
        double n = fAbs(val);
        double sum = 0.0;

        for (int i = 1; i <= n; i++) {
          sum += (Math.cos(theta(t) - t * Math.log(i))) / Math.sqrt(i);
        }
        sum = 2.0 * sum;

        double remainder;
        double frac = val - n; 
        int k = 0;
        double R = 0.0;

        // Necessary to individually calculate each remainder term by using
        // Taylor Series co-efficients. These coefficients are defined below.
        while (k <= r) {
            R = R + C(k, 2.0*frac-1.0) * Math.pow(t / twopi, 
                    ((double) k) * -0.5);
            k++;
        }

        remainder = Math.pow(-1, (int)n-1) * Math.pow(t / twopi, -0.25) * R;
        return sum + remainder;
    }

    /**
     * C terms for the Riemann-Siegel formula. See 
         * https://web.viu.ca/pughg/thesis.d/masters.thesis.pdf for details.
         * Calculates the Taylor Series coefficients for C0, C1, C2, C3, 
         * and C4. 
     * @param n - the number of coefficient terms to use.
         * @param z - referenced per the Taylor series calculations.
     * @return the Taylor series approximation of the remainder terms.
    */
    public static double C (int n, double z) {
        if (n==0) 
            return(.38268343236508977173 * Math.pow(z, 0.0) 
            +.43724046807752044936 * Math.pow(z, 2.0) 
            +.13237657548034352332 * Math.pow(z, 4.0) 
            -.01360502604767418865 * Math.pow(z, 6.0) 
            -.01356762197010358089 * Math.pow(z, 8.0) 
            -.00162372532314446528 * Math.pow(z,10.0) 
            +.00029705353733379691 * Math.pow(z,12.0) 
            +.00007943300879521470 * Math.pow(z,14.0) 
            +.00000046556124614505 * Math.pow(z,16.0) 
            -.00000143272516309551 * Math.pow(z,18.0) 
            -.00000010354847112313 * Math.pow(z,20.0) 
            +.00000001235792708386 * Math.pow(z,22.0) 
            +.00000000178810838580 * Math.pow(z,24.0) 
            -.00000000003391414390 * Math.pow(z,26.0) 
            -.00000000001632663390 * Math.pow(z,28.0) 
            -.00000000000037851093 * Math.pow(z,30.0) 
            +.00000000000009327423 * Math.pow(z,32.0) 
            +.00000000000000522184 * Math.pow(z,34.0) 
            -.00000000000000033507 * Math.pow(z,36.0) 
            -.00000000000000003412 * Math.pow(z,38.0)
            +.00000000000000000058 * Math.pow(z,40.0) 
            +.00000000000000000015 * Math.pow(z,42.0)); 
        else if (n==1) 
            return(-.02682510262837534703 * Math.pow(z, 1.0) 
            +.01378477342635185305 * Math.pow(z, 3.0) 
            +.03849125048223508223 * Math.pow(z, 5.0) 
            +.00987106629906207647 * Math.pow(z, 7.0) 
            -.00331075976085840433 * Math.pow(z, 9.0) 
            -.00146478085779541508 * Math.pow(z,11.0) 
            -.00001320794062487696 * Math.pow(z,13.0) 
            +.00005922748701847141 * Math.pow(z,15.0) 
            +.00000598024258537345 * Math.pow(z,17.0) 
            -.00000096413224561698 * Math.pow(z,19.0) 
            -.00000018334733722714 * Math.pow(z,21.0) 
            +.00000000446708756272 * Math.pow(z,23.0) 
            +.00000000270963508218 * Math.pow(z,25.0) 
            +.00000000007785288654 * Math.pow(z,27.0)
            -.00000000002343762601 * Math.pow(z,29.0) 
            -.00000000000158301728 * Math.pow(z,31.0) 
            +.00000000000012119942 * Math.pow(z,33.0) 
            +.00000000000001458378 * Math.pow(z,35.0) 
            -.00000000000000028786 * Math.pow(z,37.0) 
            -.00000000000000008663 * Math.pow(z,39.0) 
            -.00000000000000000084 * Math.pow(z,41.0) 
            +.00000000000000000036 * Math.pow(z,43.0) 
            +.00000000000000000001 * Math.pow(z,45.0)); 
      else if (n==2) 
            return(+.00518854283029316849 * Math.pow(z, 0.0) 
            +.00030946583880634746 * Math.pow(z, 2.0) 
            -.01133594107822937338 * Math.pow(z, 4.0) 
            +.00223304574195814477 * Math.pow(z, 6.0) 
            +.00519663740886233021 * Math.pow(z, 8.0) 
            +.00034399144076208337 * Math.pow(z,10.0) 
            -.00059106484274705828 * Math.pow(z,12.0) 
            -.00010229972547935857 * Math.pow(z,14.0) 
            +.00002088839221699276 * Math.pow(z,16.0) 
            +.00000592766549309654 * Math.pow(z,18.0) 
            -.00000016423838362436 * Math.pow(z,20.0) 
            -.00000015161199700941 * Math.pow(z,22.0) 
            -.00000000590780369821 * Math.pow(z,24.0) 
            +.00000000209115148595 * Math.pow(z,26.0) 
            +.00000000017815649583 * Math.pow(z,28.0) 
            -.00000000001616407246 * Math.pow(z,30.0) 
            -.00000000000238069625 * Math.pow(z,32.0) 
            +.00000000000005398265 * Math.pow(z,34.0) 
            +.00000000000001975014 * Math.pow(z,36.0) 
            +.00000000000000023333 * Math.pow(z,38.0) 
            -.00000000000000011188 * Math.pow(z,40.0) 
            -.00000000000000000416 * Math.pow(z,42.0) 
            +.00000000000000000044 * Math.pow(z,44.0) 
            +.00000000000000000003 * Math.pow(z,46.0)); 
      else if (n==3) 
            return(-.00133971609071945690 * Math.pow(z, 1.0) 
            +.00374421513637939370 * Math.pow(z, 3.0) 
            -.00133031789193214681 * Math.pow(z, 5.0) 
            -.00226546607654717871 * Math.pow(z, 7.0) 
            +.00095484999985067304 * Math.pow(z, 9.0) 
            +.00060100384589636039 * Math.pow(z,11.0) 
            -.00010128858286776622 * Math.pow(z,13.0) 
            -.00006865733449299826 * Math.pow(z,15.0) 
            +.00000059853667915386 * Math.pow(z,17.0) 
            +.00000333165985123995 * Math.pow(z,19.0)
            +.00000021919289102435 * Math.pow(z,21.0) 
            -.00000007890884245681 * Math.pow(z,23.0) 
            -.00000000941468508130 * Math.pow(z,25.0) 
            +.00000000095701162109 * Math.pow(z,27.0) 
            +.00000000018763137453 * Math.pow(z,29.0) 
            -.00000000000443783768 * Math.pow(z,31.0) 
            -.00000000000224267385 * Math.pow(z,33.0) 
            -.00000000000003627687 * Math.pow(z,35.0) 
            +.00000000000001763981 * Math.pow(z,37.0) 
            +.00000000000000079608 * Math.pow(z,39.0) 
            -.00000000000000009420 * Math.pow(z,41.0) 
            -.00000000000000000713 * Math.pow(z,43.0) 
            +.00000000000000000033 * Math.pow(z,45.0) 
            +.00000000000000000004 * Math.pow(z,47.0)); 
      else 
            return(+.00046483389361763382 * Math.pow(z, 0.0) 
            -.00100566073653404708 * Math.pow(z, 2.0) 
            +.00024044856573725793 * Math.pow(z, 4.0) 
            +.00102830861497023219 * Math.pow(z, 6.0) 
            -.00076578610717556442 * Math.pow(z, 8.0) 
            -.00020365286803084818 * Math.pow(z,10.0) 
            +.00023212290491068728 * Math.pow(z,12.0) 
            +.00003260214424386520 * Math.pow(z,14.0) 
            -.00002557906251794953 * Math.pow(z,16.0) 
            -.00000410746443891574 * Math.pow(z,18.0) 
            +.00000117811136403713 * Math.pow(z,20.0) 
            +.00000024456561422485 * Math.pow(z,22.0) 
            -.00000002391582476734 * Math.pow(z,24.0) 
            -.00000000750521420704 * Math.pow(z,26.0) 
            +.00000000013312279416 * Math.pow(z,28.0) 
            +.00000000013440626754 * Math.pow(z,30.0) 
            +.00000000000351377004 * Math.pow(z,32.0) 
            -.00000000000151915445 * Math.pow(z,34.0) 
            -.00000000000008915418 * Math.pow(z,36.0) 
            +.00000000000001119589 * Math.pow(z,38.0) 
            +.00000000000000105160 * Math.pow(z,40.0) 
            -.00000000000000005179 * Math.pow(z,42.0) 
            -.00000000000000000807 * Math.pow(z,44.0) 
            +.00000000000000000011 * Math.pow(z,46.0) 
            +.00000000000000000004 * Math.pow(z,48.0));
    }     
}

Java- Commons.Math BrentSolver returns NoBracketingException。这是一个错误吗？

Java- Commons.Math BrentSolver returns NoBracketingException. Is this a bug?

java

apache-commons

apache-commons-math