数值积分出现符号错误与精度问题 如何重写被积函数满足精度要求

我需要对如下函数进行积分：

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

我的问题是：如何重写被积函数才能达到要求的计算精度？

解决方案 #

核心问题定位 #

1. 灾难性抵消：原代码中t1 = num - den1 - den2是两个接近的大数相减，大z下arg值很高，对数项的主项会完全抵消，仅剩浮点误差，甚至出现负数结果。
2. 误差限设置不合理：原代码设置的1e-5相对误差限对1e-43量级的结果没有意义，积分器会提前终止计算。
3. 被积函数计算方式冗余：t2部分的计算存在溢出风险，没有利用渐近特性优化。

重写后的被积函数 #

#include <math.h>

double integrand__W(double x, double z) {
    const double arg = z * z / (4.0 * x);
    double t1;

    // 分区间计算log项，避免抵消
    if (arg > 30.0) { // arg很大时用渐近展开，直接计算小量
        const double exp_neg_arg = exp(-arg);
        const double exp_neg_x = exp(-x);
        t1 = exp_neg_arg * (1 + exp_neg_x) + 0.5 * exp_neg_arg * exp_neg_arg * (1 - exp(-2 * x));
    } else if (arg < 1e-3) { // arg很小时展开低阶项
        const double inv_arg = 4.0 * x / (z * z);
        t1 = log(inv_arg * (1 + exp(-x) * inv_arg)) - arg / 2.0 + arg * arg / 12.0;
    } else { // 中间区间直接计算，无抵消风险
        const double num_log = log1p(exp(-arg - x)) + arg;
        const double den_log = log(expm1(arg));
        t1 = num_log - den_log;
    }

    // 重写t2，用双曲余弦优化，避免溢出
    double t2;
    if (x > 30.0) {
        t2 = exp(-x);
    } else {
        const double cosh_half_x = cosh(0.5 * x);
        t2 = 0.25 / (cosh_half_x * cosh_half_x);
    }

    return t1 * t2;
}

积分调用优化 #

将原相对误差限改为绝对误差限，适配小量级结果的计算需求：

hcubature(1, inf_W, &params, 1, xl, xu, maxEval, 1e-50, 0, ERROR_INDIVIDUAL, &W, &W_ERR);

如果z取值继续增大，可将double替换为long double，或接入MPFR多精度计算库进一步提升精度。

内容的提问来源于stack exchange，提问作者surrutiaquir