title: Java- Looking for advice on calculating the min/max for a function or the derivative in step intervals
body: <p>Looking for advice on a math problem that has turned into a Java nightmare. I scanned the web and could not find a solution. I have looked at similar programs and unfortunately couldn't find help.</p>
  
  <p>Summary of the problem: I am looking to implement a method inside Java that will find either the min or max of the Riemann-Siegel Z(t) function (I have already created the code to calculate Z(t)) or the value of its derivative. To show what I am trying to do do, the graph of Z(t) from 0 &lt; t &lt; 100 looks like this.</p>
  
  <p><a href="https://i.stack.imgur.com/jO0WG.gif" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/jO0WG.gif" alt="enter image description here"></a></p>
  
  <p>Directly reviewing the function in <a href="http://www.wolframalpha.com/input/?i=Plot%20Z%28t%29%20from%20t%20%3D%200%20to%20t%20%3D%20100" rel="nofollow noreferrer">Wolfram Alpha</a> or in <a href="http://mathworld.wolfram.com/Riemann-SiegelFormula.html" rel="nofollow noreferrer">here</a> makes the "Java nightmare" I am having look overly complicated. The issue that I am describing is not super complicated, although it may be due to my inexperience in numerical analysis. The general outline of what I am looking to do is</p>
  
  <ol>
  <li><p>Write a method inside Java to calculate all the places where the derivative of this function is zero (in the graph above, the function has about 30 or so values between 0 &lt; t &lt; 100).</p></li>
  <li><p>Inside the method, define a step interval to evaluate the function through a lower bound and an upper bound.</p></li>
  <li><p>One of the following three methods - Calculate the max/min in one method, calculate the max/min in two methods, or calculate the values where the derivative is zero.</p></li>
  <li><p>Add this to my existing program (I have made a test program to make the problem easier. The test program looks at cos(x))</p></li>
  </ol>
  
  <p>I scanned the internet and found <a href="http://www.dsalgo.com/2013/02/FindMinMax.php.html" rel="nofollow noreferrer">this</a>. I found a lot of other different approaches but none of these seem to work. All of the solutions provided appear to calculate only one maximum/minimum/derivative inside a step interval. In order to make use of the new method, the program would need to find all the values where the derivative is zero or when the function has either a maximum or minimum. As an example, cos(x) has about 16 zeroes between 0 &lt; x &lt; 50 (the new method would find all of these values).</p>
  
  <p>To make this easier, I created a test program that can be analyzed against the cos(x) function. </p>
  
  <pre><code>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 &lt; 0.0)
              return -1;
          else if (x &gt; 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 &lt;= 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 &lt;= 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) &lt; 0.001) {
              System.out.println("The x value is " + x + ". The value of the "
                      + "derivative is " + dy);
          }
      }
  </code></pre>
  
  <p>The purpose of the test program is to check whether or not the method for <code>public static void findDerivative</code> is correct. It sort of works, although it returns two values for the approximation of the derivative. The graph of cos(x) is shown below.</p>
  
  <p><a href="https://i.stack.imgur.com/cZIo8.gif" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/cZIo8.gif" alt="enter image description here"></a> </p>
  
  <p>The values outputted by the program are </p>
  
  <pre><code>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
  </code></pre>
  
  <p>It gets close, although it needs to compute the derivative twice because of the Math.abs(dy) &lt; 0.001 inside the findDerivative method. The following ways around this have all been unsuccessful.</p>
  
  <ol>
  <li><p>A recommendation was asked to calculate the derivative through Newton's method. I don't know any way of applying Newton's because I don't know the derivative for Z(t).</p></li>
  <li><p>All of the programs I found online and in other websites directly calculate only "one" minimum or maximum in an interval from [a, b]. In the graph above and in the graph for the Z(t) function, I am looking for all the minimums and maximums (or, alternatively, when the function is zero). Calculating one minimum or maximum between an interval of [0, 100] doesn't help, I would need to have a method that will calculate all of them.</p></li>
  <li><p>I originally underestimated the difficulty of doing this.</p></li>
  </ol>
  
  <p>Does anyone have a suggestion? What could I do that would do this with the cos(x) test program? If I got this working I could go and figure out the Z(t) program myself. I have spent a lot of time thinking about this and have lost some sleep. I haven't been able to think of a way around this on my own.</p>
  
  <p>Here is what I am using to calculate the Z(t) function for general values (it is not necessary to understand the program below to work around these difficulties).</p>
  
  <pre><code>/**************************************************************************
  **
  **    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 &lt; 0.0)
              return -1;
          else if (x &gt; 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 &lt;= 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 &lt;= upperBound ; x += step) {
              if (y &gt; (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 &lt;= upperBound; i+=step) {
          double currEval = f.f(i);
          if (currEval &lt; 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 &lt;= 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 &lt;= 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));
      }     
  }
  </code></pre>