Java- Apache Math.Commons 库中的 DerivativeStructure

Java- DerivativeStructure inside the Apache Math.Commons library

我在实现 DerivativeStructure class provided by the Apache Commons Math Library 的更复杂形式时遇到问题。我写了两个程序,一个是试用程序,然后是一个真正的程序。我遇到的问题涉及我的真实程序。我认为可以通过界面解决此问题。

这个post比较长,请看我在post末尾提出的问题(我在post末尾添加了额外的信息来查明问题]).

第一个程序很简单,为 sin(x) 函数求根。

import java.util.TreeSet;
import org.apache.commons.math3.analysis.differentiation.DerivativeStructure;
import org.apache.commons.math3.analysis.differentiation.UnivariateDifferentiableFunction;
import org.apache.commons.math3.analysis.solvers.*;
import org.apache.commons.math3.exception.DimensionMismatchException;

public class Test5 {

    public static void main(String args[]) {
        NewtonRaphsonSolver test = new NewtonRaphsonSolver(1E-10);

        UnivariateDifferentiableFunction f = new UnivariateDifferentiableFunction() {

            public double value(double x) {
                return Math.sin(x);
            }

            public DerivativeStructure value(DerivativeStructure t) throws
                    DimensionMismatchException {
                return t.sin();
            }
        };

        double EPSILON = 1e-6;
        TreeSet<Double> set = new TreeSet<>();
        for (int i = 1; i <= 5000; i++) {
            set.add(test.solve(1000, f, i, i + EPSILON));
        }
        for (Double s : set) {
            if (s > 0) {
                System.out.println(s);
            }
        }

    }
}

在这个程序中,DerivativeStructureclass是通过

形成的
UnivariateDifferentiableFunction f = new UnivariateDifferentiableFunction() {

            public double value(double x) {
                return Math.sin(x);
            }

            public DerivativeStructure value(DerivativeStructure t) throws
                    DimensionMismatchException {
                return t.sin();
            }
        };

这只是 "straightforward" 因为我可以申请 ` return t.sin();'在 public DerivativeStructure 方法中。

第二个程序无法直接引用 t.sin(x) 或任何 t.function() 值,因为它是程序稍后开发的自定义函数。我尝试编写第二个程序涉及

import java.util.TreeSet;
import org.apache.commons.math3.analysis.UnivariateFunction;
import org.apache.commons.math3.analysis.differentiation.DerivativeStructure;
import org.apache.commons.math3.analysis.differentiation.
        UnivariateDifferentiableFunction;
import org.apache.commons.math3.analysis.solvers.*;
import org.apache.commons.math3.exception.DimensionMismatchException;

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

    // Main method
    public static void SiegelMain() {

        NewtonRaphsonSolver test = new NewtonRaphsonSolver(1E-10);
        UnivariateDifferentiableFunction f = new 
            UnivariateDifferentiableFunction() {

            public double value(double x) {
                return RiemennZ(x, 4);
            }

            UnivariateFunction func = (double x) -> RiemennZ(x, 4);

            public DerivativeStructure value(DerivativeStructure t) throws DimensionMismatchException {
                //I have no idea what to write here
                return new DerivativeStructure(1, 1, 0, func.value(x)));
            }

        };
        System.out.println("Zeroes inside the critical line for " +
                "Zeta(1/2 + it). The t values are referenced below.");
        System.out.println();
    double EPSILON = 1e-6;
        TreeSet<Double> set = new TreeSet<>();
        for (int i = 1; i <= 5000; i++) {
            set.add(test.solve(1000, f, i, i+1));
        }
        for (Double s : set) {
            if(s > 0)
            System.out.println(s);
        }
    }

    /**
     * Needed as a reference for the interpolation function.
    */
    public static interface Function {
    public double f(double x);
    }

    /**
     * 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 interpolation.
     * @param f - the function
         * @param lowerBound - the lower bound of integration.
         * @param upperBound - the upper bound of integration.
         * @param step - the step for dx in [a:b]
     * @return the roots of the specified function.
    */
    public static void findRoots(Function f, double lowerBound,
                  double upperBound, double step) {
    double x = lowerBound, next_x = x;
    double y = f.f(x), next_y = y;
    int s = sign(y), next_s = s;

    for (x = lowerBound; x <= upperBound ; x += step) {
        s = sign(y = f.f(x));
        if (s == 0) {
        System.out.println(x);
        } else if (s != next_s) {
        double dx = x - next_x;
        double dy = y - next_y;
        double cx = x - dx * (y / dy);
        System.out.println(cx);
        }
        next_x = x; next_y = y; next_s = 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));
    }     
}

第二个程序的问题是由于

UnivariateDifferentiableFunction f = new 
            UnivariateDifferentiableFunction() {

            public double value(double x) {
                return RiemennZ(x, 4);
            }

            UnivariateFunction func = (double x) -> RiemennZ(x, 4);

            public DerivativeStructure value(DerivativeStructure t) throws DimensionMismatchException {
                //I have no idea what to write here
                return new DerivativeStructure(1, 1, 0, func.value(x)));
            }

        };

具体涉及实施

public DerivativeStructure value(DerivativeStructure t) throws DimensionMismatchException {
                //I have no idea what to write here
                return new DerivativeStructure(1, 1, 0, func.value(x)));
            }

在 Apache Commons 数学库中,this 页末尾列出了以下文档。

用户可以通过多种方式创建 UnivariateDifferentiableFunction 接口的实现。第一种方法是简单地直接使用 DerivativeStructure 中的适当方法来计算加法、减法、正弦、余弦……这通常非常简单,无需记住微分规则:用户代码仅表示函数本身,微分将在引擎盖下自动计算。第二种方法是编写 classical UnivariateFunction 并将其传递给 UnivariateFunctionDifferentiator 接口的现有实现,以检索同一函数的微分版本。第一种方法更适合用户已经控制所有底层代码的小功能。第二种方法更适合使用 DerivativeStructure API 编写起来很麻烦的大型函数,或者用户无法控制完整底层代码的函数(例如调用外部库的函数)。

有没有办法通过编写然后通过接口实现 DerivativeStructure 来解决这个问题?

我不太清楚,RiemennZ 是什么,但它似乎是某种函数族。在这种情况下,您应该以为第二个参数的每个值产生一个新的 UnivariateDifferentiableFunction 的方式实现它:

public static UnivariateDifferentiableFunction RiemennZ(int r) {
  return new UnivariateDifferentiableFunction() {
    double value(double t) {
       /* YOUR CODE FROM ABOVE */
    }
    public DerivativeStructure value(DerivativeStructure t) {
       /* YOUR CODE FROM ABOVE BUT AS A DS */
    }
}

这称为 Higher order function,在其他语言中很常见。您必须仅为您实际想要区分的参数提供 DerivativeStructure 实现。