Java- Looking for advice on calculating the min/max for a function or the derivative in step intervals
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.
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 < t < 100 looks like this.
Directly reviewing the function in Wolfram Alpha or in here 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
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 < t < 100).
Inside the method, define a step interval to evaluate the function through a lower bound and an upper bound.
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.
Add this to my existing program (I have made a test program to make the problem easier. The test program looks at cos(x))
I scanned the internet and found this. 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 < x < 50 (the new method would find all of these values).
To make this easier, I created a test program that can be analyzed against the cos(x) function.
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);
}
}
The purpose of the test program is to check whether or not the method for public static void findDerivative 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.
The values outputted by the program are
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
It gets close, although it needs to compute the derivative twice because of the Math.abs(dy) < 0.001 inside the findDerivative method. The following ways around this have all been unsuccessful.
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).
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.
I originally underestimated the difficulty of doing this.
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.
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).
/**************************************************************************
**
** 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)
+.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));
}
}
Solution
It looks like you're trying to do numerical optimization. The Apache Commons Math library has several implementations for both optimization and root-finding. Even if you ultimately have to write your own implementation, it might be helpful to first prototype your solution using the algorithms available in the library to find one that works before you implement it yourself.
- Is there a java.time equivalent of DateTime.dayOfYear().isLeap()?
- WordNet Java API
- Why use a ReentrantLock if one can use synchronized(this)?
- How to format a ZonedDateTime to a String?
- How to short circuit when using a parallel stream with filter and findFirst
- Stream distinct with nested list counts and matching nested list item counts
- Timezone changes are not being persisted in Spring-Boot application
- how to convert decimal to hexadecimal in java
- Why does pattern matching with switch on InetAddress fail with 'does not cover all possible input values'?
- Jackson to parse bean<t> with list of T beans
- Catch a generic exception in Java?
- How to get the whole text of StAX XMLEvent object?
- The www-authenticate header breaks JettyClientHttpConnector
- Why should I use @NaturalId?
- Can't reproduce virtual thread pinning (jdk 21), yet mysql not as parallel as should be
- how to move slider towards left and right using selenium java
- Format of Date in the JavaFX.TableView
- Launch Jetty 9 Programmatically
- Spring JdbcTemplate fail to convert String to Enum (that extend specific interface), since default converter is called even if removed
- GET Request Returning Empty Object in SpringBoot
- Trying to draw box around checkbox and display text using pdfbox
- Extracting date, hour and minute from Instant.now() generated date
- Tail Call Optimisation in Java
- Spring Boot Angular18 Uploading Photo
- What is a IOException, and how do I fix it?
- What does %5B and %5D in POST requests stand for?
- Why couldn't Integer type convert to Double type automatically in Java8 stream API?
- Spring Boot default H2 jdbc connection (and H2 console)
- No content to map due to end-of-input jackson parser
- How to properly add HTTP header in WebService Client (Spring Web Services)