数值积分中的符号错误，可能存在精度问题

Question

我需要集成以下功能：

其中 z > 0。问题是被积函数对于大 z 来说非常小，并且在积分时需要高精度。到目前为止，我已经将被积函数写为

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;
}

对于数值积分，我使用 Cubature，一个用于自适应多维积分的简单 C 程序包：

//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;

其中通过变量的变化可以在半无限区间上积分：

通过解析，我知道积分是非负的，所以积分本身应该是非负的。但是，由于缺乏准确性，我得到了一些不正确的结果：

z: 100 | -3.97632e-17 , 1.24182e-16

在Mathematica中，以高精度工作，我可以得到想要的结果：

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 *)

我的问题：有什么方法可以写出我的被积函数以达到要求的精度吗？谢谢。

Answer 1

关于数学我不能说太多（我与数学有 love/hate 关系）但是更高的精度可以通过 long double 和标准数学库中的相关数学函数。

但 long double 并不一定意味着更高的精度，取决于您的编译器和系统架构，它可能只是双精度或 80 位扩展精度或更高。

更多信息：

Answer 2

有些积分方案只使用正权重，如果被积函数的评估函数值都是非负的，就会产生非负的积分值。其他一些积分方案允许负权重，从而可能获得更高的积分精度。 Cubature 可能使用其中之一。

对于 z=100，您的实际积分值非常接近 0，这也是您得到的结果，因此积分方案确实没有任何问题。如果您绝对需要非负性，一种选择是将负结果简单地设置为 0。

Answer 3

在向另一个 community 询问相同的问题后，我得到了两个似乎有效的建议：

避免减法抵消

先稍微操作一下积分：

$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)}$

然后将被积函数重写为

double integrand__W(double x, double z){
    double arg = z*z/(4.0*x);
    
    double t1 = log1p((exp(-x)+1)/expm1(arg));
    
    double num = exp(x);
    double den = exp(x)+1;
    double t2 = isinf(den) ? exp(-x) : num/(den*den);
    
    return t1*t2;
}

Exp-Sinh 正交的使用

此集成方案由Boost库提供：

#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";
}

数值积分中的符号错误，可能存在精度问题

Wrong sign in numerical integration, possible precision issue

c

numerical-integration

避免减法抵消

Exp-Sinh 正交的使用