2015-08-29 7 views
4

Ищите совет по математической проблеме, превратившейся в кошмар Java. Я просмотрел Интернет и не смог найти решение. Я посмотрел похожие программы и, к сожалению, не смог найти помощь.Java- Ищет совет по вычислению min/max для функции или производной с интервалом шага

Резюме проблемы: Я ищу для реализации метода внутри Java, который найдет либо min, либо max функции Riemann-Siegel Z (t) (я уже создал код для вычисления Z (t)) или значение его производной. Чтобы показать, что я пытаюсь сделать, график Z ​​(t) от 0 < t < 100 выглядит следующим образом.

enter image description here

Непосредственно рассматривает функцию в Wolfram Alpha или here делает «Java кошмар» Я имею взгляд слишком сложный. Проблема, которую я описываю, не слишком сложна, хотя это может быть из-за моей неопытности в численном анализе. Общая схема того, что я ищу, чтобы сделать

  1. Напишите метод внутри Java, чтобы вычислить все места, где производная этой функции равна нулю (в приведенном выше графике, функция имеет около 30 или около того значения между 0 < t < 100).

  2. Внутри метода задайте интервал шагов для оценки функции через нижнюю границу и верхнюю границу.

  3. Один из следующих трех методов. Рассчитайте максимальное/минимальное значение в одном методе, вычислите max/min двумя способами или вычислите значения, где производная равна нулю.

  4. Добавить это мой существующую программу (я сделал тестовую программу, чтобы сделать проблему проще. Тестовая программа смотрит на соз (х))

Я просмотрел в Интернете и нашел this. Я нашел много других подходов, но ни один из них, похоже, не работает. Все предоставленные решения, по-видимому, вычисляют только один максимум/минимальный/производный в интервале шагов. Чтобы использовать новый метод, программе нужно было бы найти все значения, где производная равна нулю или когда функция имеет максимум или минимум. Например, cos (x) имеет около 16 нулей между 0 < x < 50 (новый метод найдет все эти значения).

Чтобы сделать это проще, я создал тестовую программу, которая может быть проанализирована относительно функции 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). При приближении х к 0 существует бесконечно много экстремумов. Для произвольных функций невозможно узнать, когда вы «делаете» поиск экстремумов на некотором интервале. Функция может выглядеть гладкой при выборке в одном масштабе, но в более мелком масштабе может быть что-то скрывающее. Вот почему большинство математических библиотек просто возвращают один макс или мин и называют это днем. – Sam

ответ

1

Похоже, вы пытаетесь выполнить численную оптимизацию. Библиотека Apache Commons Math имеет несколько реализаций для optimization и root-finding. Даже если вы в конечном счете должны написать свою собственную реализацию, было бы полезно сначала прототипировать ваше решение, используя алгоритмы, доступные в библиотеке, чтобы найти тот, который работает до его реализации.