gsl Fehler im inifinit Integrationsintervall. Schlechtes integriertes Verhalten gefunden. Wie man es repariert?

Question

Ich erhalte die folgende Fehlermeldung nach dem Versuch, die eine numerische Integration auf einem infinte Intervall [0, inf) unter Verwendung von GSL in C.

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

Hier wird die Funktion Ich Integration $

zu tun

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

Und der Code für das Integrationsverfahren

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

die fucntion den folgenden Aspekt hat:

Welche, für meine Augen, hat ein sehr gutes Verhalten. Anstatt also eine unendliche Integration durchgeführt wird, fahre ich die Integration bis zu einer Obergrenze, die ich für rezonable, wie in:

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

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

ein Ergebnis zu erzielen, die für unendliche Integration mit dem Ergebnis von Mathematica stimmt

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

Aber das GSL unendliche integrale keps ist "Nan". Irgendwelche Ideen? Ich danke im Voraus für die Hilfe.

Answer 1

Ich denke, das Problem hier ist, dass C im Gegensatz zu Mathematica keine willkürliche Genauigkeit in der Berechnung verwendet. Dann, zu einem Zeitpunkt, an dem Exp [Ep] berechnet wird, läuft die numerische Berechnung über. Jetzt

, GSL verwendet die Transformation x = (1-t)/t, auf Abstand zur Karte (0,1]. also für t < < 0 posible ist, da das Verhalten Ihrer Funktion nan Ergebnisse zu erhalten neigt dazu, Unbestimmtheiten (0/0 oder inf/inf, usw.) für extreme Werte Vielleicht, wenn Sie die Bedingungen

Exp schreiben [(Ep (x) - \ Mu)/T]./{1 + Exp [ (Ep (x) - \ Mu)/T]}^2

Verwendung A/B = Exp [Ln. A - Ln B], könnten Sie ein besseren numerisches Verhalten bekommen

ich versuchen, Wenn ich gute Ergebnisse habe, werde ich es dir sagen.

Die Lösung

Wie ich schon sagte, um die Probleme kümmern entstehen mit unbestimmten Formen annehmen müssen. So können die problematischen Bedingungen schreiben die logarithmische Version mit:

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

und mit dieser Hauptfunktion

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

Sie erhalten die Antwort: -5.065288e-03. Ich bin neugierig ...Dies ist, wie ich die Funktion in Mathematica definieren

So vergleicht die Antworten:

• GSL -5.065288e-03
• Mathematica -0,005065287633739702

Answer 2

Als @ Yonatan zuleta ochoa weist richtig darauf hin, das Problem ist in exp(t)/pow(exp(t)+1,2). exp(t) kann ein ieee754 DBL_MAX für Werte von t so niedrig wie nextafter(log(DBL_MAX), INFINITY) überlaufen, was ~7.09783e2 ist.

Wenn exp(t) == INFINITY,

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

Yonatan vorgeschlagene Lösung ist Logarithmen zu verwenden, was getan werden kann, wie folgt:

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 

Dies ist ein ganz vernünftiger Ansatz, NAN bis zu sehr hohen Werten zu vermeiden von t. Jedoch in Ihrem Code kann t == (Ep ± mu)/TINFINITY sein, wenn abs(T) < 1.0 für Werte von x der Nähe von DBL_MAX, auch wenn x ist nicht Unendlichkeit. In diesem Fall wird die Subtraktion t - 2*log1p(exp(t)) zu ∞ - ∞, was wiederum NAN ist.

Ein anderer Ansatz wird durch Dividieren sowohl Nenner und Zähler durch exp(x) (die nicht null für jedes finite x ist) exp(x)/pow(exp(x)+1,2) mit 1.0/(pow(exp(x)+1,2)*pow(exp(x), -1)) zu ersetzen. Dies vereinfacht sich zu 1.0/(exp(x)+exp(-x)+2.0). Hier

ist eine Implementierung der Funktion NAN für Werte von x bis zu vermeiden und einschließlich DBL_MAX:

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

Diese Implementierung auch hypot(x,M) verwendet, anstatt sqrt(x*x, M*M) und vermeidet x*x Berechnung durch die Reihenfolge der Multiplikationen Neuanordnung/Abteilungen zu Gruppe x/Ep zusammen. Da hypot(x,M)abs(x) für abs(x) >> abs(M) sein wird, nähert sich der Begriff x/Ep1.0 für große x.