(Note: I found a partial solution. I have pasted the test results at the end)
I am looking to numerically integrate an approximation of the zeta function. This is done through what is known as the Abel-Plana formula. The Abel-Plana formula can be used to numerically evaluate Zeta(s) for s = a + i*b.
I first wrote a program to calculate Zeta(s) for s in real. It seems to work,
/**************************************************************************
**
** Abel-Plana Formula for the Zeta Function
**
**************************************************************************
** Axion004
** 08/16/2015
**
** This program computes the value for Zeta(s) using a definite integral
** approximation through the Abel-Plana formula. The Abel-Plana formula
** can be shown to approximate the value for Zeta(s) through a definite
** integral. The integral approximation is handled through the Composite
** Simpson's Rule known as Adaptive Quadrature.
**************************************************************************/
import java.util.*;
import java.math.*;
public class AbelMain2 {
public static void main(String[] args) {
AbelMain();
}
public static void AbelMain() {
double s = 0;
double start, stop, totalTime;
Scanner scan = new Scanner(System.in);
System.out.print("Enter the value of s inside the Riemann Zeta " +
"Function: ");
try {
s = scan.nextDouble();
}
catch (Exception e) {
System.out.println("Please enter a valid number for s.");
}
start = System.currentTimeMillis();
System.out.println("The value for Zeta(s) is " + AbelPlana(s));
stop = System.currentTimeMillis();
totalTime = (double) (stop-start) / 1000.0;
System.out.println("Total time taken is " + totalTime + " seconds.");
}
// The definite integral for Zeta(s) in the Abel Plana formula.
// Numerator = Sin(s * arctan(t))
// Denominator = (1 + t^2)^(s/2) * (e^(pi*t) + 1)
public static double function(double x, double s) {
double num = Math.sin(s * Math.atan(x));
double den = Math.pow((1.0 + Math.pow(x, 2.0)), s / 2.0) *
(Math.pow(Math.E, Math.PI * x) + 1.0);
return num / den;
}
// Adaptive quadrature - See http://www.mathworks.com/moler/quad.pdf
// Used to approximate values for definite integrals
// Parameters - a is the start of the integral, b is the end of the
// integral, s is the double value used to evaulate Zeta(s).
public static double adaptiveQuad(double a, double b, double s) {
double EPSILON = 1E-6;
double step = b - a;
double c = (a + b) / 2.0;
double d = (a + c) / 2.0;
double e = (b + c) / 2.0;
double S1 = step / 6.0 * (function(a, s) + 4*function(c, s) +
function(b, s));
double S2 = step / 12.0 * (function(a, s) + 4*function(d, s) +
2*function(c, s) + 4*function(e, s) + function(b, s));
if (Math.abs(S2 - S1) <= EPSILON)
return S2 + (S2 - S1) / 15.0;
else
return adaptiveQuad(a, c, s) + adaptiveQuad(c, b, s);
}
// Returns an approximate sum of Zeta(s) through an integral aproximation
// by Adaptive Quadrature
public static double AbelPlana(double s) {
double C1 = Math.pow(2.0, s - 1.0) / (s - 1.0);
double C2 = Math.pow(2.0, s);
return C1 - C2 *adaptiveQuad(0, 10, s);
}
}
I then tried to extend this program to s in complex. It appears that the program doesn't work because I am trying to numerically evaluate a complex function through quadrature. The quadrature method that I wrote is only applicable to real-value functions as shown in Adaptive Quadrature. The work around I used is rather strange although it seems to work.
Here is what I wrote for the numerical approximation. I had to use a helper class in order to calculate complex numbers in Java.
/**************************************************************************
**
** Abel-Plana Formula for the Zeta Function
**
**************************************************************************
** Axion004
** 08/16/2015
**
** This program computes the value for Zeta(s) using a definite integral
** approximation through the Abel-Plana formula. The Abel-Plana formula
** can be shown to approximate the value for Zeta(s) through a definite
** integral. The integral approximation is handled through the Composite
** Simpson's Rule known as Adaptive Quadrature.
**************************************************************************/
import java.util.*;
import java.math.*;
public class AbelMain extends Complex {
public static void main(String[] args) {
AbelMain();
}
public static void AbelMain() {
double re = 0;
double im = 0;
double start, stop, totalTime;
Scanner scan = new Scanner(System.in);
System.out.println("Calculation of the Riemann Zeta " +
"Function in the form Zeta(s) = a + ib.");
System.out.println();
System.out.print("Enter the value of [a] inside the Riemann Zeta " +
"Function: ");
try {
re = scan.nextDouble();
}
catch (Exception e) {
System.out.println("Please enter a valid number for a.");
}
System.out.print("Enter the value of [b] inside the Riemann Zeta " +
"Function: ");
try {
im = scan.nextDouble();
}
catch (Exception e) {
System.out.println("Please enter a valid number for b.");
}
start = System.currentTimeMillis();
Complex z = new Complex(re, im);
if ( re == 1 && im == 0)
System.out.println("Zeta(s) is undefined for Zeta(1 + 0*i).");
else
System.out.println("The value for Zeta(s) is " + AbelPlana(z));
stop = System.currentTimeMillis();
totalTime = (double) (stop-start) / 1000.0;
System.out.println("Total time taken is " + totalTime + " seconds.");
}
// The definite integral for Zeta(s) in the Abel Plana formula.
// Numerator = Sin(s * arctan(t))
// Denominator = (1 + t^2)^(s/2) * (e^(pi*t) + 1)
public static Complex f(double t, Complex z) {
Complex num = (z.multiply(Math.atan(t))).sin();
Complex D1 = new Complex(1 + t*t, 0).pow(z.divide(2.0));
double D2 = Math.pow(Math.E, Math.PI * t) + 1.0;
Complex den = D1.multiply(D2);
return num.divide(den);
}
// Adaptive quadrature - See http://www.mathworks.com/moler/quad.pdf
// Used to approximate values for definite integrals
// Parameters - a is the start of the integral, b is the end of the
// integral, s is the double value used to evaulate Zeta(s).
public static Complex adaptiveQuad(double a, double b, Complex s) {
double EPSILON = 1E-6;
double step = b - a;
double c = (a + b) / 2.0;
double d = (a + c) / 2.0;
double e = (b + c) / 2.0;
Complex S1 = (f(a, s).add(f(c, s).multiply(4)).add(f(b, s))).
multiply(step / 6.0);
Complex S2 = (f(a, s).add(f(d, s).multiply(4)).add(f(c, s).multiply(2))
.add(f(e, s).multiply(4)).add(f(b, s))).multiply(step / 12.0);
Complex result = (S2.minus(S1)).divide(15.0);
if(S2.minus(S1).mod() <= EPSILON)
return S2.add(result);
else
return adaptiveQuad(a, c, s).add(adaptiveQuad(c, b, s));
}
// Returns an approximate sum of Zeta(s) through an integral aproximation
// by Adaptive Quadrature
public static Complex AbelPlana(Complex z) {
Complex two = new Complex(2.0, 0.0);
Complex C1 = two.pow(z.minus(1.0)).divide(z.minus(1.0));
Complex C2 = two.pow(z);
Complex mult = C2.multiply(adaptiveQuad(0, 10, z));
if ( z.re < 0 && z.re % 2 == 0 && z.im == 0)
return new Complex(0.0, 0.0);
else
return C1.minus(mult);
}
public AbelMain(double re, double im) {
super(re, im);
}
}
The second class for Complex numbers.
/**************************************************************************
**
** Complex Numbers
**
**************************************************************************
** Axion004
** 08/20/2015
**
** This class is necessary as a helper class for the calculation of
** imaginary numbers. The calculation of Zeta(s) inside AbelMain is in
** the form of z = a + i*b.
**************************************************************************/
public class Complex extends Object{
public double re;
public double im;
/**
Constructor for the complex number z = a + i*b
@param re Real part
@param im Imaginary part
*/
public Complex (double re, double im) {
this.re = re;
this.im = im;
}
/**
Real part of the Complex number
@return Re[z] where z = a + i*b.
*/
public double real() {
return re;
}
/**
Imaginary part of the Complex number
@return Im[z] where z = a + i*b.
*/
public double imag() {
return im;
}
/**
Complex conjugate of the Complex number
in which the conjugate of z is z-bar.
@return z-bar where z = a + i*b and z-bar = a - i*b
*/
public Complex conjugate() {
return new Complex(re, -im);
}
/**
Addition of two Complex numbers (z is unchanged).
<br>(a+i*b) + (c+i*d) = (a+c) + i*(b+d)
@param t is the number to add.
@return z+t where z = a+i*b and t = c+i*d
*/
public Complex add(Complex t) {
return new Complex(re + t.real(), im + t.imag());
}
/**
Addition of Complex number and a double.
@param d is the number to add.
@return z+d where z = a+i*b and d = double
*/
public Complex add(double d){
return new Complex(this.re + d, this.im);
}
/**
Subtraction of two Complex numbers (z is unchanged).
<br>(a+i*b) - (c+i*d) = (a-c) + i*(b-d)
@param t is the number to subtract.
@return z-t where z = a+i*b and t = c+i*d
*/
public Complex minus(Complex t) {
return new Complex(re - t.real(), im - t.imag());
}
/**
Subtraction of Complex number and a double.
@param d is the number to subtract.
@return z-d where z = a+i*b and d = double
*/
public Complex minus(double d){
return new Complex(this.re - d, this.im);
}
/**
Complex multiplication (z is unchanged).
<br> (a+i*b) * (c+i*d) = (a*c)+ i(b*c) + i(a*d) - (b*d)
@param t is the number to multiply by.
@return z*t where z = a+i*b and t = c+i*d
*/
public Complex multiply(Complex t) {
return new Complex(re * t.real() - im * t.imag(), re *
t.imag() + im * t.real());
}
/**
Complex multiplication by a double.
@param d is the double to multiply by.
@return z*d where z = a+i*b and d = double
*/
public Complex multiply(double d){
return new Complex(this.re * d,this.im * d);}
/**
Modulus of a Complex number or the distance from the origin in
* the polar coordinate plane.
@return |z| where z = a + i*b.
*/
public double mod() {
if ( re != 0.0 || im != 0.0)
return Math.sqrt(re*re + im*im);
else
return 0.0;
}
/**
* Modulus of a Complex number squared
* @param z = a + i*b
* @return |z|^2 where z = a + i*b
*/
public double abs(Complex z) {
return Math.sqrt(Math.pow(re*re, 2) + Math.pow(im*im, 2));
}
/**
Division of Complex numbers (doesn't change this Complex number).
<br>(a+i*b) / (c+i*d) = (a+i*b)*(c-i*d) / (c+i*d)*(c-i*d) =
* (a*c+b*d) + i*(b*c-a*d) / (c^2-d^2)
@param t is the number to divide by
@return new Complex number z/t where z = a+i*b
*/
public Complex divide (Complex t) {
double denom = Math.pow(t.mod(), 2); // Square the modulus
return new Complex((re * t.real() + im * t.imag()) / denom,
(im * t.real() - re * t.imag()) / denom);
}
/**
Division of Complex number by a double.
@param d is the double to divide
@return new Complex number z/d where z = a+i*b
*/
public Complex divide(double d){
return new Complex(this.re / d, this.im / d);
}
/**
Exponential of a complex number (z is unchanged).
<br> e^(a+i*b) = e^a * e^(i*b) = e^a * (cos(b) + i*sin(b))
@return exp(z) where z = a+i*b
*/
public Complex exp () {
return new Complex(Math.exp(re) * Math.cos(im), Math.exp(re) *
Math.sin(im));
}
/**
The Argument of a Complex number or the angle in radians
with respect to polar coordinates.
<br> Tan(theta) = b / a, theta = Arctan(b / a)
<br> a is the real part on the horizontal axis
<br> b is the imaginary part of the vertical axis
@return arg(z) where z = a+i*b.
*/
public double arg() {
return Math.atan2(im, re);
}
/**
The log or principal branch of a Complex number (z is unchanged).
<br> Log(a+i*b) = ln|a+i*b| + i*Arg(z) = ln(sqrt(a^2+b^2))
* + i*Arg(z) = ln (mod(z)) + i*Arctan(b/a)
@return log(z) where z = a+i*b
*/
public Complex log() {
return new Complex(Math.log(this.mod()), this.arg());
}
/**
The square root of a Complex number (z is unchanged).
Returns the principal branch of the square root.
<br> z = e^(i*theta) = r*cos(theta) + i*r*sin(theta)
<br> r = sqrt(a^2+b^2)
<br> cos(theta) = a / r, sin(theta) = b / r
<br> By De Moivre's Theorem, sqrt(z) = sqrt(a+i*b) =
* e^(i*theta / 2) = r(cos(theta/2) + i*sin(theta/2))
@return sqrt(z) where z = a+i*b
*/
public Complex sqrt() {
double r = this.mod();
double halfTheta = this.arg() / 2;
return new Complex(Math.sqrt(r) * Math.cos(halfTheta), Math.sqrt(r) *
Math.sin(halfTheta));
}
/**
The real cosh function for Complex numbers.
<br> cosh(theta) = (e^(theta) + e^(-theta)) / 2
@return cosh(theta)
*/
private double cosh(double theta) {
return (Math.exp(theta) + Math.exp(-theta)) / 2;
}
/**
The real sinh function for Complex numbers.
<br> sinh(theta) = (e^(theta) - e^(-theta)) / 2
@return sinh(theta)
*/
private double sinh(double theta) {
return (Math.exp(theta) - Math.exp(-theta)) / 2;
}
/**
The sin function for the Complex number (z is unchanged).
<br> sin(a+i*b) = cosh(b)*sin(a) + i*(sinh(b)*cos(a))
@return sin(z) where z = a+i*b
*/
public Complex sin() {
return new Complex(cosh(im) * Math.sin(re), sinh(im)*
Math.cos(re));
}
/**
The cos function for the Complex number (z is unchanged).
<br> cos(a +i*b) = cosh(b)*cos(a) + i*(-sinh(b)*sin(a))
@return cos(z) where z = a+i*b
*/
public Complex cos() {
return new Complex(cosh(im) * Math.cos(re), -sinh(im) *
Math.sin(re));
}
/**
The hyperbolic sin of the Complex number (z is unchanged).
<br> sinh(a+i*b) = sinh(a)*cos(b) + i*(cosh(a)*sin(b))
@return sinh(z) where z = a+i*b
*/
public Complex sinh() {
return new Complex(sinh(re) * Math.cos(im), cosh(re) *
Math.sin(im));
}
/**
The hyperbolic cosine of the Complex number (z is unchanged).
<br> cosh(a+i*b) = cosh(a)*cos(b) + i*(sinh(a)*sin(b))
@return cosh(z) where z = a+i*b
*/
public Complex cosh() {
return new Complex(cosh(re) *Math.cos(im), sinh(re) *
Math.sin(im));
}
/**
The tan of the Complex number (z is unchanged).
<br> tan (a+i*b) = sin(a+i*b) / cos(a+i*b)
@return tan(z) where z = a+i*b
*/
public Complex tan() {
return (this.sin()).divide(this.cos());
}
/**
The arctan of the Complex number (z is unchanged).
<br> tan^(-1)(a+i*b) = 1/2 i*(log(1-i*(a+b*i))-log(1+i*(a+b*i))) =
<br> -1/2 i*(log(i*a - b+1)-log(-i*a + b+1))
@return arctan(z) where z = a+i*b
*/
public Complex atan(){
Complex ima = new Complex(0.0,-1.0); //multiply by negative i
Complex num = new Complex(this.re,this.im-1.0);
Complex den = new Complex(-this.re,-this.im-1.0);
Complex two = new Complex(2.0, 0.0); // divide by 2
return ima.multiply(num.divide(den).log()).divide(two);}
/**
* The Math.pow equivalent of two Complex numbers.
* @param z - the complex base in the form z = a + i*b
* @return z^y where z = a + i*b and y = c + i*d
*/
public Complex pow(Complex z){
Complex a = z.multiply(this.log());
return a.exp();
}
/**
* The Math.pow equivalent of a Complex number to the power
* of a double.
* @param d - the double to be taken as the power.
* @return z^d where z = a + i*b and d = double
*/
public Complex pow(double d){
Complex a=(this.log()).multiply(d);
return a.exp();
}
/**
Override the .toString() method to generate complex numbers, the
* string representation is now a literal Complex number.
@return a+i*b, a-i*b, a, or i*b as desired.
*/
public String toString() {
if (re != 0.0 && im > 0.0) {
return re + " + " + im +"i";
}
if (re !=0.0 && im < 0.0) {
return re + " - "+ (-im) + "i";
}
if (im == 0.0) {
return String.valueOf(re);
}
if (re == 0.0) {
return im + "i";
}
return re + " + i*" + im;
}
public static void main(String[] args) {
Complex s = new Complex(2.0, 3.0);
Complex y = new Complex(3.0, 4.0);
Complex d = new Complex(4.0, 2.0);
Complex x = new Complex(1.0, 0.0);
System.out.println(s.atan());
System.out.println(s.divide(d));
System.out.println(s.im);
System.out.println(s.pow(y));
System.out.println(s.pow(2.0));
}
}
I did some preliminary testing and it seems that both the public static Complex f and public static Complex AbelPlana methods work fine. I checked these calculations against Mathematica. The method for
public static Complex adaptiveQuad(double a, double b, Complex z)
also works (for low values).
Calculation of the Riemann Zeta Function in the form Zeta(s) = a + ib.
Enter the value of [a] inside the Riemann Zeta Function: -12
Enter the value of [b] inside the Riemann Zeta Function: 1.2
The value for Zeta(s) is 0.0900630360334386 + 0.08241454022912213*i
Total time taken is 0.014 seconds.
Enter the value of [a] inside the Riemann Zeta Function: 0.3
Enter the value of [b] inside the Riemann Zeta Function: -2.1
The value for Zeta(s) is 0.4003421605328948 + 0.2570024810252463*i
Total time taken is 0.005 seconds.
Enter the value of [a] inside the Riemann Zeta Function: 1
Enter the value of [b] inside the Riemann Zeta Function: 2
The value for Zeta(s) is 0.598165521081844 - 0.35185471257583556*i
Total time taken is 0.005 seconds.
Enter the value of [a] inside the Riemann Zeta Function: 0.5
Enter the value of [b] inside the Riemann Zeta Function: 10
The value for Zeta(s) is 1.5448963436469043 - 0.1153378574814412*i
Total time taken is 0.014 seconds.
Something is clearly wrong with larger values.
Enter the value of [a] inside the Riemann Zeta Function: 24
Enter the value of [b] inside the Riemann Zeta Function: -8
The value for Zeta(s) is 164561.70457820524 + 302628.94685036544*i
Total time taken is 0.003 seconds.
I scanned the internet and couldn't find a method for numerical integration of complex functions. The integral for the approximation is shown in here. It doesn't seem like I can use this relationship in the program (it needs to calculate the Zeta(z) for z = a+i*b anyhow).
I searched and found this, and this. Do I need to reference the Apache Commons Math library to numerically integrate a complex function? I would prefer to write a method to do this myself. My current knowledge doesn't seem to be sufficient, all of the numerical integration methods that I know directly approximate real-valued functions.
I will review why my method fails for large values...
