Numerical integration giving incorrect sign

Question

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,
  Method -> "LocalAdaptive"]

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.

Answer 1

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>
#include <boost/math/quadrature/exp_sinh.hpp>

using boost::math::quadrature::exp_sinh;
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/math/quadrature/exp_sinh.hpp>
#include <boost/multiprecision/float128.hpp>
using boost::multiprecision::float128;
using boost::math::quadrature::exp_sinh;
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";
}
```

Answer 2

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?

Answer 3

(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.