Apache Math.Commons库中的JavaDerivativeStructure [英] Java- DerivativeStructure inside the Apache Math.Commons library

查看:85
本文介绍了Apache Math.Commons库中的JavaDerivativeStructure的处理方法,对大家解决问题具有一定的参考价值,需要的朋友们下面随着小编来一起学习吧!

问题描述

我在实现 Apache Commons Math Library提供的.html"rel =" nofollow> DerivativeStructure类.我编写了两个程序,一个是试用程序,然后是一个真实程序.我遇到的问题与我的真实程序有关.我认为可以通过界面解决该问题.

I am having trouble implementing a more complicated form of the DerivativeStructure class provided by the Apache Commons Math Library. I wrote two programs, a trial program and then a real program. The problem I have encountered deals with my real program. I think that it is possible to work around the issue through an interface.

这篇文章相当长,请查看我在文章末尾提出的问题(我在文章末尾添加了更多信息以查明问题).

This post is rather long, please see the problem I have presented at the end of the post (I added extra information to pinpoint the issue at the end of the post).

第一个程序很简单,它为sin(x)函数找到了根.

The first program is straightforward and finds roots for the sin(x) function.

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

    }
}

在此程序中,DerivativeStructure类通过以下方式形成

Inside this program, the DerivativeStructure class is formed through

UnivariateDifferentiableFunction f = new UnivariateDifferentiableFunction() {

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

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

这只是直截了当",因为我可以申请`返回t.sin();'在公共DerivativeStructure方法中.

This is only "straightforward" because I am able to apply ` return t.sin();' inside the public DerivativeStructure method.

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

The second program is not able to directly reference t.sin(x) or any t.function() value because it is a custom function that is developed later in the program. My attempt at writing the second program involves

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

第二个程序的问题是由于

The problem with the second program is due to

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 Math库中,

Inside the Apache Commons Math library, the following documentation is listed at the end of this page.

用户可以通过多种方式创建UnivariateDifferentiableFunction接口的实现.第一种方法是直接使用DerivativeStructure的适当方法直接编写它,以计算加法,减法,正弦,余弦...这通常是很简单的,不需要记住区分规则:用户代码仅表示功能本身,差异将在引擎盖下自动计算.第二种方法是编写经典的UnivariateFunction并将其传递给UnivariateFunctionDifferentiator接口的现有实现,以检索同一函数的差异版本.第一种方法更适合于用户已经控制了所有基础代码的小型功能.第二种方法更适合使用DerivativeStructure API编写笨重的大型函数,或者用户无法控制全部底层代码的函数(例如,调用外部库的函数).

是否有某种方法可以解决此问题,方法是编写并通过接口实现DerivativeStructure?

Is there some way of working around this issue through writing and then implementing the DerivativeStructure by an interface?

推荐答案

我不完全知道 RiemennZ 是什么,但它似乎是某种函数家族.在这种情况下,您应该以对第二个参数的每个值产生一个新的 UnivariateDifferentiableFunction 的方式来实现它:

I do not know exactly, what RiemennZ is, but it seems to be some kind of function family. In that case, you should implement it in a way that yields a new UnivariateDifferentiableFunction for each value of the second parameter:

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

这称为高阶函数,在其他语言中也很常见.您只需为您实际想要区分的参数提供 DerivativeStructure 实现.

This is called a Higher order function, and quite common in other languages. You have to provide the DerivativeStructure implementation only for the argument you actually want to differentiate.

这篇关于Apache Math.Commons库中的JavaDerivativeStructure的文章就介绍到这了,希望我们推荐的答案对大家有所帮助,也希望大家多多支持IT屋!

查看全文
登录 关闭
扫码关注1秒登录
发送“验证码”获取 | 15天全站免登陆