I need to integrate the following function:
where z > 0. The problem is that the integrand is very small for large z and high precision is required in the integration. So far, I have written the integrand as
double integrand__W(double x, double z){
double arg = z*z/(4.0*x);
double num = exp(arg x) 1;
double den1 = expm1(arg);
double den2 = exp(x);
num = isinf(num) ? arg x : log(num);
den1 = isinf(den1) ? arg : log(den1);
den2 = x; //log(exp(x))=x
double t1 = num-den1-den2;
num = exp(x);
double den = exp(x) 1;
double t2 = isinf(den) ? exp(-x) : num/(den*den);
return t1*t2;
}
For numerical integration, I'm using 
Analytically, I know that the integrated is non-negative, so the integral itself should be non-negative. However, I'm getting some incorrect results due to a lack of accuracy:
z: 100 | -3.97632e-17 , 1.24182e-16
In Mathematica, working with high precision, I can get the desired result:
w[x_, z_] := E^x/(E^x 1)^2 Log[(E^(z^2/(4 x)) E^-x)/(E^(z^2/(4 x)) - 1)]
W[z_?NumericQ] := NIntegrate[w[x, z], {x, 0, ∞},
WorkingPrecision -> 40,
Method -> "LocalAdaptive"]
W[100]
(* 4.679853458969239635780655689865016458810*10^-43 *)
My question: Is there any way to write my integrand such that I can reach the required precision? Thanks.
CodePudding user response：
There are integration schemes which only use positive weights, resulting in nonnegative integral values if the evaluated function values of the integrand are all nonnegative. Some other integration schemes permit negative weights, resulting in a possibly higher accuracy for integration. Cubature probably uses one of those.
Your actual integral value is very close to 0 for z=100, and that's what you're getting, too, so there's really nothing wrong with the integration scheme. If you absolutely need nonnegativity, one option is to simply set the negative results to 0.
CodePudding user response：
I can't say much about the mathematics (i have a love/hate relationship with math) but higher precision can be achieved via long double and the related math functions in the standard math library.
But long double does not necessarily mean higher precision, dependend on your compiler and system architecture it could be simply a double or 80 bit extended precision or more.
More info: