我在实现 DerivativeStructure class 提供的更复杂形式的 Apache Commons Math Library 时遇到问题。我写了两个程序，一个是试用程序，然后是一个真正的程序。我遇到的问题与我的实际程序有关。我认为可以通过界面解决这个问题。
这篇文章相当长，请参阅我在文章末尾提出的问题(我在文章末尾添加了额外的信息来查明问题)。
第一个程序很简单，它找到了 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);
}
}
}
}
在这个程序中，DerivativeStructure类是通过形成的
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();'在公共(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 Math 库中，this 页面末尾列出了以下文档。
用户可以通过多种方式创建 UnivariateDifferentiableFunction 接口(interface)的实现。第一种方法是简单地使用 DerivativeStructure 中的适当方法直接编写它来计算加法、减法、正弦、余弦...这通常非常简单，无需记住微分规则:用户代码仅代表函数本身，差异将在后台自动计算。第二种方法是编写经典的 UnivariateFunction 并将其传递给 UnivariateFunctionDifferentiator 接口(interface)的现有实现，以检索同一函数的微分版本。第一种方法更适合用户已经控制所有底层代码的小功能。第二种方法更适合使用 DerivativeStructure API 编写起来很麻烦的大型函数，或者用户无法控制完整底层代码的函数(例如调用外部库的函数)。
是否有某种方法可以通过编写然后通过接口(interface)实现 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 实现。
关于Apache Math.Commons 库中的 Java-DerivativeStructure，我们在Stack Overflow上找到一个类似的问题： https://stackoverflow.com/questions/32511458/