Java - численное интегрирование комплексной функции - дзета-функция, формула Абеля-Планы
(Примечание: я нашел частичное решение. В конце я вставил результаты теста)
Я хочу численно интегрировать приближение дзета-функции. Это делается с помощью так называемой формулы Абеля-Планы. Формула Абеля-Планы может быть использована для численной оценки дзетов для s = a + i*b.
Сначала я написал программу для вычисления Zeta(s) для s в реальном времени. Похоже на работу,
/**************************************************************************
**
** 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);
}
}
Затем я попытался расширить эту программу до сложных. Похоже, что программа не работает, потому что я пытаюсь численно оценить сложную функцию с помощью квадратуры. Квадратурный метод, который я написал, применим только к функциям с действительными значениями, как показано в Адаптивной квадратуре. Работа вокруг, которую я использовал, довольно странная, хотя кажется, что она работает.
Вот что я написал для числового приближения. Мне пришлось использовать вспомогательный класс для вычисления комплексных чисел в 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);
}
}
Второй класс для комплексных чисел.
/**************************************************************************
**
** 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));
}
}
Я провел предварительное тестирование, и кажется, что как общедоступные статические методы Complex f, так и общедоступные статические методы Complex AbelPlana работают нормально. Я проверил эти расчеты с Mathematica. Метод для
public static Complex adaptiveQuad(double a, double b, Complex z)
также работает (для низких значений).
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.
Я просканировал Интернет и не смог найти метод численного интегрирования сложных функций. Интеграл для приближения показан здесь. Похоже, я не могу использовать это соотношение в программе (ему все равно нужно вычислить дзета (z) для z = a+i*b).
Я искал и нашел это, и это. Нужно ли ссылаться на библиотеку Apache Commons Math для численной интеграции сложной функции? Я бы предпочел написать метод, чтобы сделать это сам. Мое текущее знание не кажется достаточным, все методы численного интегрирования, которые я знаю, непосредственно приближают к вещественным функциям.
Я расскажу, почему мой метод не работает для больших значений...