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

Posted 2023-03-24

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

【英文标题】：Wrong sign in numerical integration, possible precision issue

【发布时间】：2021-12-24 23:35:57

【问题描述】：

我需要集成以下功能：

在哪里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】：

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

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

【参考方案2】：

在向不同的community 提出相同的问题后，我得到了两条似乎可行的建议：

避免减法抵消

先对积分进行一点操作：

然后将被积函数改写为

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

【参考方案3】：

我不能对数学说太多（我对数学有爱/恨的关系），但通过long double 和标准数学库中的相关数学函数可以实现更高的精度 .

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

更多信息：

https://en.wikipedia.org/wiki/Long_double https://en.wikipedia.org/wiki/Extended_precision

【讨论】：

提高精度并不能解决集成方案中的问题。