2015-08-29 51 views
4

尋找已經變成Java惡夢的數學問題的建議。我掃描了網頁並找不到解決方案。我看過類似的計劃,很遺憾找不到幫助。 (我已經創建了計算Z(t)的代碼)的最小值或最大值。我想要在Java中實現一個方法,它可以找到Riemann-Siegel Z(t)函數的最小值或最大值(我已經創建了計算Z(t)的代碼)或其衍生物的價值。爲了顯示我想要做的事情,從0 < t < 100的Z(t)的圖看起來像這樣。Java-尋找關於計算函數的最小/最大值或步長間隔的導數的建議

enter image description here

直接審查功能Wolfram Alphahere使「Java的噩夢」我有看過於複雜。我所描述的問題並不是非常複雜,儘管這可能是由於我在數值分析方面缺乏經驗。什麼我希望做的大致輪廓是

  1. 寫Java裏面的方法來計算所有的地方,這個函數的導數爲零(在上面的圖中,功能大約有30個左右的值在0 < t < 100)之間。

  2. 在該方法中,定義一個步驟間隔以通過下界和上界來評估函數。

  3. 以下三種方法之一 - 計算一種方法中的最大/最小值,計算兩種方法中的最大/最小值,或計算導數爲零時的值。

  4. 這添加到我的現有程序(我已經做了一個測試程序,以使問題更容易。該測試程序會在cos(x))

我掃描上網,發現this。我發現了很多其他不同的方法,但這些方法似乎都不起作用。所提供的所有解決方案似乎只能在步驟間隔內計算一個最大值/最小值/導數。爲了使用新方法,程序需要找到所有導數爲零的值,或者當函數具有最大值或最小值時。例如,cos(x)在0 < x < 50之間有大約16個零(新方法可以找到所有這些值)。

爲了使這更容易,我創建了一個測試程序,可以根據cos(x)函數進行分析。

import java.math.*; 

public class Test { 
    public static void main(String[] args){ 
     Function cos = new Function() 
     { 
     public double f(double x) { 
     return Math.cos(x); 
     } 
    }; 


     //findRoots(cos, 1, 1000, 0.001); 
     findDerivative(cos, 1, 100, 0.001); 
    } 

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

    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 the specified function passed in with a lower bound, 
    // upper bound, and step size. 
    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; 
    } 
    } 

    public static void findDerivative(Function f, double lowerBound, 
      double upperBound, double step) { 

    for (double x = lowerBound; x <= upperBound; x += step) { 
     double fxstep = f.f(x); 
     double fx = fxstep; 
     fxstep = f.f(x+step); 
     double dy = (fxstep - fx)/step; 
     if (Math.abs(dy) < 0.001) { 
      System.out.println("The x value is " + x + ". The value of the " 
        + "derivative is " + dy); 
     } 
    } 

測試程序的目的是檢查public static void findDerivative該方法是否是正確的。它有點類似,但它返回了導數近似值的兩個值。 cos(x)的圖如下所示。

enter image description here

由程序輸出的值是

The x value is 3.140999999999764. The value of the derivative is -9.265358602572604E-5 
The x value is 3.141999999999764. The value of the derivative is 9.073462475805982E-4 
The x value is 6.282000000000432. The value of the derivative is 6.853070969592423E-4 
The x value is 6.283000000000432. The value of the derivative is -3.1469280259432963E-4 
The x value is 9.424000000000216. The value of the derivative is -2.7796075396935294E-4 
The x value is 9.425000000000216. The value of the derivative is 7.220391380347024E-4 
The x value is 12.564999999998475. The value of the derivative is 8.706142144987439E-4 
The x value is 12.565999999998475. The value of the derivative is -1.2938563354047972E-4 
The x value is 15.706999999996734. The value of the derivative is -4.632679163618647E-4 
The x value is 15.707999999996733. The value of the derivative is 5.36731999623008E-4 
The x value is 18.849000000000053. The value of the derivative is 5.592153640154862E-5 
The x value is 18.850000000000055. The value of the derivative is -9.440782817726756E-4 
The x value is 21.990000000003892. The value of the derivative is -6.485750521090239E-4 
The x value is 21.991000000003893. The value of the derivative is 3.514248534397524E-4 
The x value is 25.132000000007732. The value of the derivative is 2.4122869812792658E-4 
The x value is 25.133000000007733. The value of the derivative is -7.587711848833223E-4 
The x value is 28.27300000001157. The value of the derivative is -8.338821652076334E-4 
The x value is 28.274000000011572. The value of the derivative is 1.6611769582119962E-4 
The x value is 31.41500000001541. The value of the derivative is 4.2653585174967645E-4 
The x value is 31.416000000015412. The value of the derivative is -5.734640621257725E-4 
The x value is 34.55700000001016. The value of the derivative is -1.9189476674341677E-5 
The x value is 34.55800000001016. The value of the derivative is 9.808103242914257E-4 
The x value is 37.69800000000284. The value of the derivative is 6.118430110335638E-4 
The x value is 37.69900000000284. The value of the derivative is -3.881568994001938E-4 
The x value is 40.83999999999552. The value of the derivative is -2.0449666182642545E-4 
The x value is 40.84099999999552. The value of the derivative is 7.955032111928162E-4 
The x value is 43.9809999999882. The value of the derivative is 7.971501513326373E-4 
The x value is 43.9819999999882. The value of the derivative is -2.028497212425151E-4 
The x value is 47.12299999998088. The value of the derivative is -3.8980383987308187E-4 
The x value is 47.123999999980875. The value of the derivative is 6.10196070671698E-4 
The x value is 50.26399999997356. The value of the derivative is 9.824572642092022E-4 
The x value is 50.264999999973554. The value of the derivative is -1.754253620145363E-5 
The x value is 53.405999999966234. The value of the derivative is -5.75111004597062E-4 
The x value is 53.40699999996623. The value of the derivative is 4.2488890927838696E-4 
The x value is 56.54799999995891. The value of the derivative is 1.6776464961676396E-4 
The x value is 56.54899999995891. The value of the derivative is -8.322352119671805E-4 
The x value is 59.68899999995159. The value of the derivative is -7.604181495590723E-4 
The x value is 59.68999999995159. The value of the derivative is 2.39581733230132E-4 
The x value is 62.83099999994427. The value of the derivative is 3.530718295507995E-4 
The x value is 62.831999999944266. The value of the derivative is -6.469280763310437E-4 
The x value is 65.97199999995095. The value of the derivative is -9.457252543310091E-4 
The x value is 65.97299999995096. The value of the derivative is 5.4274563066059045E-5 
The x value is 69.11399999996596. The value of the derivative is 5.383789610791112E-4 
The x value is 69.11499999996596. The value of the derivative is -4.616209549057615E-4 
The x value is 72.25599999998096. The value of the derivative is -1.3103257845425986E-4 
The x value is 72.25699999998096. The value of the derivative is 8.689672701400752E-4 

它靠攏,儘管它需要計算導數的兩倍,因爲Math.abs(DY)< 0.001 findDerivative方法內的。以下方法都是不成功的。

  1. 推薦被要求通過牛頓法計算導數。我不知道應用牛頓的任何方法,因爲我不知道Z(t)的導數。

  2. 我在網上和其他網站上找到的所有程序都直接計算從[a,b]開始的最小或最大間隔。在上圖和Z(t)函數圖中,我正在尋找所有的最小值和最大值(或者,當函數爲零時)。計算間隔[0,100]之間的一個最小值或最大值不會有幫助,我需要一個方法來計算所有這些值。

  3. 我原本低估了這樣做的難度。

有沒有人有建議?我可以做什麼,這將與cos(x)測試程序做到這一點?如果我有這個工作,我可以自己去弄清楚Z(t)程序。我花了很多時間思考這個問題,並失去了一些睡眠。我一直無法想象解決這個問題的方法。

下面是我用來計算一般值的Z(t)函數(沒有必要了解下面的程序來解決這些困難)。

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

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

    // Main method 
    public static void SiegelMain() { 
     Function RiemennSiegelZ = new Function() 
     { 
     public double f(double x) { 
     return RiemennZ(x, 4); 
     } 
    }; 
     System.out.println("Zeroes inside the critical line for " + 
       "Zeta(1/2 + it). The t values are referenced below."); 
     System.out.println(); 
     // Uncomment to find non-trivial zeroes for Zeta(1/2 + it) 
    findRoots(RiemennSiegelZ, 1, 40000, 0.001); 
     //findMax(RiemennSiegelZ, 1, 400, 0.001); 
    } 

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

    /** 
    * Calculates the local maximum from a provided lower and upper bound. 
    * @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 local maximum for the function. 
    */ 
    public static void findMax(Function f, double lowerBound, 
        double upperBound, double step) { 
    double x = lowerBound, next_x = x + step; 
    double y = f.f(x), next_y = y + step; 

    for (x = lowerBound; x <= upperBound ; x += step) { 
      if (y > (next_y)) { 
     System.out.println(y); 
     } 
     next_x = x; next_y = y; 
    } 
    } 

    /** 
    * Calculates the local minimum from a provided lower and upper bound. 
    * @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 local minimum for the function. 
    */ 
    public static double findMin(Function f, double lowerBound, double 
      upperBound, double step) { 
    double minValue = f.f(lowerBound); 

    for (double i=lowerBound; i <= upperBound; i+=step) { 
     double currEval = f.f(i); 
     if (currEval < minValue) { 
      minValue = currEval; 
     } 
    } 

     return minValue; 
    } 

    /** 
    * 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) 
      +.000000* 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)); 
    }  
} 
+0

你有沒有試過一些Java數學庫? – nashuald

+0

我知道數學的公共圖書館。我期望自己實施這個。我認爲自己解決會更有益處。 – Axion004

+0

作爲一個思想實驗,考慮你的程序在提供像sin(1/x)這樣的函數時應該做些什麼。當x接近0時,存在無限多的極值。對於任意函數,無法知道在某個時間間隔內何時「完成」尋找極值。當以一個比例採樣時,函數可能看起來很平滑,但是在更精細的範圍內可以隱藏任何東西。這就是爲什麼大多數數學庫只返回一個最大值或最小值並稱之爲一天。 – Sam

回答

1

它看起來像你試圖做數值優化。 Apache Commons Math庫有optimizationroot-finding的幾個實現。即使您最終必須編寫自己的實現,首先使用庫中可用的算法對您的解決方案進行原型化,以便在您自己實現之前找到可用的解決方案可能會有所幫助。