gsl無限長積分區間中的錯誤。發現不好的被積函數行爲。如何解決它？

Question

我試圖做的C.

gsl: qags.c:553: ERROR: bad integrand behavior found in the integration interval 
Default GSL error handler invoked. 
Command terminated by signal 6 

這裏使用GSL一個infinte區間[0，INF）的數值積分後，得到以下錯誤消息，我的集成功能 $

double dI2dmu(double x, void * parametros){ 
    double *p,Ep,mu,M,T; 
    p=(double *) parametros; 

    M=p[0]; 
    T=p[1]; 
    mu=p[2]; 

    Ep=sqrt(x*x+M*M); 

    double fplus= -((exp((Ep - mu)/T)/(pow(1 + exp((Ep - mu)/T),2)*T) - exp((Ep + \ 
mu)/T)/(pow(1 + exp((Ep + mu)/T),2)*T))*pow(x,2))/(2.*Ep*pow(PI,2)); 
    return fplus; 
} 

而對於整合過程

params[0]=0.007683; //M 
params[1]=0.284000;// T 
params[2]=0.1; //mu 

    gsl_function dI2mu_u; 
    dI2mu_u.function = &dI2dmu; 
    dI2mu_u.params = &params; 
    gsl_integration_qagiu (&dI2mu_u, 0, 0, 1e-7, 100000, 
      w, &resultTest2, &error1Test2); 

的溫控功能具有以下方面的代碼：

對我而言，它的行爲非常好。所以，與其進行無限融合，我進行集成高達，我認爲rezonable，就像一個上限：

gsl_function G; 
G.function = &dI2dmu; 
G.params = &params; 

gsl_integration_qags (&G, 0, 1e2*A, 0, 1e-7, 100000, 
        w, &result1, &error1); 

獲取與數學的無限融合的結果一致的結果

result definite up to 10*A  = 0.005065263943958745 
result up to infinity    = nan 
Mathematica result up to infinity = 0.005065260000000000 

但GSL無限積分keps爲「nan」。有任何想法嗎？我預先感謝您的幫助。

Answer 1

我認爲這裏的問題是，不像Mathematica，C在計算中不使用任意精度。然後，在計算Exp [Ep]時的某個點，數值計算溢出。

現在，GSL使用變換X =（1-T）/ T，到映射到間隔（0,1]。 所以，對於t < < 0是更多鈔票獲得因爲你的函數的行爲楠導致往往會出現極端值的不確定性（0/0或inf/inf等） 也許如果你寫出條款

Exp [（Ep（x） - \ Mu）/ T]/{1 + Exp [ （EP（X） - \畝）/ T]}^2

使用A/B = EXP [LN A - 。LN B]，你可以得到一個更好的數值行爲

我會努力如果和我有很好的結果，那麼我會告訴你。

解決方案

正如我之前說的，你必須要注意不確定的形式出現的問題。所以，讓我們寫出來的使用對數版本的問題條款：

double dIdmu(double x, void * parametros){ 
     double *p,Ep,mu,M,T; 
     p=(double *) parametros; 

     M=p[0]; 
     T=p[1]; 
     mu=p[2]; 

     Ep=sqrt(x*x+M*M); 

    double fplus= - (exp((Ep - mu)/T -2.0*log(1.0 + exp((Ep - mu)/T))) - exp((Ep + mu)/T -2.0*log(1.0 + exp((Ep + mu)/T)))) * pow(x,2) /(2.* T * Ep*pow(M_PI,2)); 

return fplus; 
     } 

，並用此功能主要

int main() 
{ 
    double params[3]; 

    double resultTest2, error1Test2; 

    gsl_integration_workspace * w 
    = gsl_integration_workspace_alloc (10000); 

    params[0]=0.007683; //M 
    params[1]=0.284000;// T 
    params[2]=0.1; //mu 

    gsl_function dI2mu_u; 
    dI2mu_u.function = &dIdmu; 
    dI2mu_u.params = &params; 
    gsl_integration_qagiu (&dI2mu_u, 0.0, 1e-7, 1e-7, 10000, w, &resultTest2, &error1Test2); 


    printf("%e\n", resultTest2); 
    gsl_integration_workspace_free (w); 

    return 0; 
} 

你得到的答案是： -5.065288e-03。 我很好奇......這是我如何定義功能在數學

所以比較答案：

• GSL -5.065288e-03
• 數學-0.005065287633739702

Answer 2

正如@ yonatan zuleta ochoa指出正確，問題在exp(t)/pow(exp(t)+1,2)。 exp(t)可以溢出ieee754​​的值爲t，低至nextafter(log(DBL_MAX), INFINITY)，即~7.09783e2。

當exp(t) == INFINITY，

exp(t)/pow(exp(t)+1,2) == ∞/pow(∞+1,2) == ∞/∞ == NAN 

Yonatan提出的解決方案如下使用對數，這是可以做到：

exp(t)/pow(exp(t)+1,2) == exp(log(exp(t)) - log(pow(exp(t)+1,2))) 
         == exp(t - 2*log(exp(t)+1)) 
         == exp(t - 2*log1p(exp(t))) //<math.h> function avoiding loss of precision for log(exp(t)+1)) if exp(t) << 1.0 

這是一個完全合理的做法，避免NAN到非常高的值t。但是，在您的代碼中，t == (Ep ± mu)/T可以是INFINITY如果abs(T) < 1.0的值爲x接近​​，即使x是而不是無限。在這種情況下，減法t - 2*log1p(exp(t))變爲∞ - ∞，這又是NAN。

一種不同的方法是通過將兩個分母和分子通過exp(x)（其不爲零對於任何有限x）取代exp(x)/pow(exp(x)+1,2)與1.0/(pow(exp(x)+1,2)*pow(exp(x), -1))。這簡化爲1.0/(exp(x)+exp(-x)+2.0)。

下面是函數避免NAN最多的x值的實施方案，並且包括​​：

static double auxfun4(double a, double b, double c, double d) 
{ 
    return 1.0/(a*b+2.0+c*d); 
} 
double dI2dmu(double x, void * parametros) 
{ 
    double *p = (double *) parametros; 
    double invT = 1.0/p[1]; 
    double Ep = hypot(x, p[0]); 
    double muexp = exp(p[2]*invT); 
    double Epexp = exp(Ep*invT); 
    double muinv = 1.0/muexp; 
    double Epinv = 1.0/Epexp; 
    double subterm = auxfun4(Epexp, muinv, Epinv, muexp); 
    subterm -= auxfun4(Epexp, muexp, Epinv, muinv); 
    double fminus = subterm*(x/Ep)*invT*(0.5/(M_PI*M_PI))*x;; 
    return -fminus; 
} 

此實現還使用hypot(x,M)，而不是sqrt(x*x, M*M)，並且避免了通過重新排列的乘法的順序計算x*x /部門將x/Ep組合在一起。由於hypot(x,M)對於abs(x) >> abs(M)將是abs(x)，所以對於大x，術語x/Ep接近1.0。