c++ - 数值积分给出了与解析表达式非常不同的结果

Question

我正在尝试计算本参考文献中定义的一些热平均积分。为了便于讨论，我们假设一个数量的平均值X如下所示：

其中M和T是参数。使用Cubature，一个用于自适应多维集成的简单 C 包，我能够实现这些积分。

最简单的情况X(k)=1具有以下分析近似值，我想用我的数值积分进行交叉检查：

其中K2是修正的贝塞尔函数。Mathematica证实了近似。这些数值积分的实现（见下文）似乎适用于虚拟示例：

./main_nNeq 30 100
0.3 | 1.77268e+06 | 1.95712e+06

但我的实际代码需要非常极端的值，其中两个值完全不同：

/main_nNeq 1e12 7.11258e17
1.40596e-06 | 4.92814e+46 | 7.19634e+53

问题：这里的潜在问题可能是什么？谢谢！

---

我的代码（无特殊原因用 C++ 编写）如下所示：

//
// COMPILING INSTRUCTIONS: 
// g++ -o main_nNeq nNeq.cpp cubature-master/hcubature.c -lgsl -lm -lgslcblas -lgmp -std=c++11
//
// MORE INFO: http://ab-initio.mit.edu/wiki/index.php/Cubature_(Multi-dimensional_integration)
//

#include <stdio.h>
#include <iostream>
#include <iomanip>
#include <fstream>
#include <string>
#include <chrono>
#include <cmath>
#include <complex>
#include "cubature-master/cubature.h"
#include <gsl/gsl_sf_dilog.h>
#include <gsl/gsl_math.h>
#include <vector>
#include <algorithm>
#include <iterator>
#include <boost/math/special_functions/bessel.hpp>


using namespace std;


#define SQR(x) ((x)*(x)) //x^2
#define CUB(x) ((x)*(x)*(x)) //x^3
#define K1(x) (boost::math::cyl_bessel_k(1.0, x)) //BesselK(1, x)
#define K2(x) (boost::math::cyl_bessel_k(2.0, x)) //BesselK(2, x)

//Momentum grid
const double log_kmin = -2;
const double log_kmax = 25;
const int Ngrid = 2000;

//Numerical constants
const double gw = 2.0;
const double PI = M_PI;
const double a_int = 0.0;

//Cosmological parameters
const double g_star = 106.75;
const double Mpl = 1.22e19; //GeV
const double aR = Mpl/2.0*sqrt(45.0/CUB(PI)/g_star); //as in Eq. (83), arXiv:1812.02651
const double Tcom = aR;
#define aa(eta) (aR*eta) //as in Eq. (82), arXiv:1812.02651


//f_F
double f_F(long double k, long double T, long double M){
    long double root = SQR(k) + SQR(M);
    root = isinf(root) ? k : sqrt(root);
    long double expo = exp(root/T);
    return isinf(expo) ? 0.0 : 1.0/(expo+1.0);
}

//n_N_eq
long double n_N_eq(double T, double M){
    return SQR(M)*T*K2(M/T);
}

//integrand
double integrand__n_N_eq(double k, double T, double M){
    return SQR(k)*f_F(k, T, M);
}

//integrator
struct fparams {
    double M;
    double T;
};

//function to be integrated
int inf_n_N_eq(unsigned ndim, const double *x, void *fdata, unsigned fdim, double *fval){
    struct fparams * fp = (struct fparams *)fdata;
    //(void)(dim); /* avoid unused parameter warnings */
    //(void)(params);
    double M = fp->M;
    double T = fp->T;
    double t = x[0]; 
    double aux = integrand__n_N_eq(a_int+t*pow(1.0-t, -1.0), T, M)*pow(1.0-t, -2.0);
    if (!isnan(aux) && !isinf(aux))
    {
        fval[0] = aux;
    }
    else
    {
        fval[0] = 0.0;
    }
    return 0;
}





int main (int argc, char* argv[])
{
    //Defining variables (M, T)
    double M = stof(argv[1]); //command line argument
    double T = stof(argv[2]); //command line argument

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

    double nNeq, nNeq_ERR;

    struct fparams params_nNeq = {M, T};
    hcubature(1, inf_n_N_eq, &params_nNeq, 1, xl, xu, 0, 0, 1e-4, ERROR_INDIVIDUAL, &nNeq, &nNeq_ERR);
    
    cout << M/T << " | " << nNeq << " | " << n_N_eq(T, M) << '\n';

    return 0;
}

标签： c++numerical-methodsnumerical-integration