# Numerical integration giving incorrect sign

For my research, I need to integrate the following function:

$$W(z)=\int_0^{\infty}dx\ w(x,z)\\ =\int_0^{\infty}dx\frac{e^x}{(e^x+1)^2}\log{\left(\frac{e^{z^2/4x+x}+1}{e^{z^2/4x+x}-e^x}\right)}\\ =\int_0^{\infty}dx\frac{e^x}{(e^x+1)^2}\log{\left(\frac{e^{z^2/4x}+e^{-x}}{e^{z^2/4x}-1}\right)}$$

where $$z > 0$$. The numerical integration of this function is part of a larger code, where the results of $$W(z)$$ are used as input for other functions. 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;
}


and its plot is shown below:

For numerical integration, I'm using Cubature, a simple C-package for adaptive multidimensional integration:

//integrator
struct fparams {
double z;
};

int inf_W(unsigned ndim, const double *x, void *fdata, unsigned fdim, double *fval){
struct fparams * fp = (struct fparams *)fdata;
double z = fp->z;
double t = x[0];
double aux = integrand__W(a_int+t*pow(1.0-t, -1.0), z)*pow(1.0-t, -2.0);
if (!isnan(aux) && !isinf(aux))
{
fval[0] = aux;
}
else
{
fval[0] = 0.0;
}
return 0;
}

//range integration 1D
size_t maxEval = 1e7;
double xl[1] = { 0 };
double xu[1] = { 1 };

double W, W_ERR;
struct fparams params = {z};
hcubature(1, inf_W, &params, 1, xl, xu, maxEval, 0, 1e-5, ERROR_INDIVIDUAL, &W, &W_ERR);
cout << "z: " << z << " | " << W << " , " << W_ERR << endl;


where the integration over the semi-infinite interval is possible by a change of variables:

$$\int_a^{\infty}f(x)dx=\int_0^1 f\left(a+\frac{t}{1-t}\right)\frac{1}{(1-t)^2}dt$$

Analytically, I know that $$w(x,z)$$ 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,

W[100]

(* 4.679853458969239635780655689865016458810*10^-43 *)


So far, I have heard two possible diagnostics regarding this issue:

1. The integrand is very small for large $$z$$ and high precision is required in the integration.
2. The integration scheme permits negative weights, sometimes resulting in possibly higher accuracy for integration, but here is causing the false negatives.

My question: Any ideas about how to fix this issue, either by rewriting my integrand such that I can reach the required precision or switching to a different integration scheme? I need to solve this integral within my code, so using Mathematica is not an option.

• It might be an idea to compute the log of this integral - you would get something that might behave better for values very close to 0 (but positive). I'm not sure if it's easy to transform to a log integral though. Another idea is to use a series expansion for large $z$, Mathematica can provide asymptotic expansions for instance.
– PC1
Nov 15, 2021 at 3:55
• @surrutiaquir Does the following help: t1=log1p((exp(-x)+1)/expm1(arg))? Nov 15, 2021 at 9:56
• I am a bit confused here.... the first formular does not match the C implementation. I believe there is an error. The error is in the expansion of the summation of the denominator of the logarithm .
– Bort
Nov 15, 2021 at 13:00
• @Bort I did a conversion from image to a MathJax formula...but I don't think I made a mistake (from the image point of view)... — surrutiaquir, can you please, double-check your original formulas, as well as my conversion? Nov 15, 2021 at 15:39
• @AntonMenshov the expression is correct, thanks! Nov 15, 2021 at 16:11

If you can use Boost, you'll get the benefit of exp-sinh quadrature, which appears to work well on this problem:

#include <iostream>
#include <cmath>

using std::exp;
using std::expm1;
using std::log;

int main() {
exp_sinh<double> integrator;
double z = 100.0;
auto f = [z](double x) {
double k1 = 1.0/(2 + exp(-x) +exp(x));
double t = z*z/(4*x);
double log_arg;
if (t > 1) {
log_arg = (1 + exp(-x)*exp(-t))/(1 - exp(-t));
} else {
log_arg = (exp(t) + exp(-x))/expm1(t);
}
return k1*log(log_arg);
};
double termination = sqrt(std::numeric_limits<double>::epsilon());
double error;
double L1;
double Q = integrator.integrate(f, termination, &error, &L1);
std::cout << "Q = " << Q << ", error estimate: " << error << "\n";
}


Output:

Q = 2.56597e-45, error estimate: 6.30636e-46


If I change double to boost::multiprecision::float128, I get

Q = 4.67986e-43, error estimate: 3.59871e-49


Edit: The diff to make it to boost::multiprecision::float128:

#include <iostream>
#include <cmath>
#include <boost/multiprecision/float128.hpp>
using boost::multiprecision::float128;
using std::exp;
using std::expm1;
using std::log;

int main() {
exp_sinh<float128> integrator;
float128 z = 100;
auto f = [z](float128 x) {
float128 k1 = 1/(2 + exp(-x) +exp(x));
float128 t = z*z/(4*x);
float128 log_arg;
if (t > 1) {
log_arg = (1 + exp(-x)*exp(-t))/(1 - exp(-t));
} else {
log_arg = (exp(t) + exp(-x))/expm1(t);
}
return k1*log(log_arg);
};
float128 termination = sqrt(std::numeric_limits<float128>::epsilon());
float128 error;
float128 L1;
float128 Q = integrator.integrate(f, termination, &error, &L1);
std::cout << "Q = " << Q << ", error estimate: " << error << "\n";
}
$$$$


This might be an accuracy problem in computing the second term, because of those large exponentials when $$x \gg 1$$.

I would first work on that term: gather $$e^x$$ out from numerator and denominator and simplify it out (leaving a summand $$e^{-x}$$ in the numerator).

However, rather than guessing it might be necessary to understand better where the accuracy problems are: can you plot the integrand $$f(x)$$ for various values of $$x$$? And compare it to the "exact" values computed by Mathematica in higher precision?

• Thanks, Federico. I have updated my post including a plot of the integrand, as well as the suggestions made by you and @njuffaabout rewriting the integrand. It seems this way works better, but I'm not sure why. Nov 15, 2021 at 18:33
• @surrutiaquir My suggested rewrite of the computation of t1` is advantageous because it avoids the issue of subtractive cancellation (sometimes called catastrophic cancellation, although I consider that a misnomer because it is not always catastrophic). Nov 15, 2021 at 21:28

(I don't have any experience with Cubature, I don't know which algorithms they implement.)

If the integration rule you use has negative weights, like Newton-Cotes rules, then even if the integrand is strictly positive, the result you get may be negative. This was an issue we would commonly encounter when using MFEM to compute norms (probably other FEM packages will have the same problem). However, in those cases, you can rigorously show that only the sign is wrong. So an easy fix for us was to take the absolute value of the result.

A solution for you might be to use a quadrature rule with positive weights, like one of many Gaussian quadrature rules.