GSL Incorrect values with numerical Integration

Question

I have an Integrand which reads

static double Integrand(double k , void * params)
{
 long double *p = (long double *)params;
  long double MassSquared = p[0];
  long double Temp = p[1];
  long double sign = p[2];

  long double EXP = std::exp(-std::sqrt(k*k+MassSquared/(Temp*Temp)));
  double res = 0;
  if(sign == -1)
  {
    res = k*k*std::log(1-exp(-std::sqrt(k*k+MassSquared/(Temp*Temp))));
  }
  else if(sign == +1)
  {
    res = k*k*std::log(1+exp(-std::sqrt(k*k+MassSquared/(Temp*Temp))));
  }
  return res;
};

This one is in a Class called VEffClass.

My Integration Routine would then be

long double VEffClass::NumIntInf(long double low, double sign, long double
MassSquared, long double Temp)
{
//long double  res = 0;
double up =
std::sqrt(std::log(C_threshold)*std::log(C_threshold)- MassSquared/(Temp*Temp));
long double StepSize = (up-low)/C_IntSteps;
long double par[3] = {MassSquared, Temp, sign};


 if(Temp == 0) return 0;

gsl_integration_workspace *w = gsl_integration_workspace_alloc(5e+5);
double result, error;
gsl_function F;
F.function = &VEffClass::Integrand;
F.params = &par;

gsl_integration_qags(&F, C_threshold, up, 1e-07, 0, 5e+5, w, &result,
&error);
gsl_integration_workspace_free(w);

return result;

}

with C_threshold = 1e-3;

If I run NumIntInf(10^-3, +1 , 100, 500) the result is 1.83038. The Result given by Maple is ~ 6*10^9 . If I use a simple Simpsons 1/3 Method it differes about 1% from the Maple result.

Can anyone point me to the Error?

Thanks for any suggestions

You need to simplify your explanation. Write down exactly what integral you want to compute and what are the parameter values used in your example. It seems you are integrating an exponential decay ( std::log(1 - exp(-..) ) - and if so you need to be careful when using an adaptive integrator. If the integral x-range is too wide in comparison to the range where the integrand is non-zero, then the integrator will think the integrand is zero! — Vivian Miranda, Jul 23 '15 at 19:53
Hello, The Integral is given by the integral over the function Integrand , theoretically from 0 to infinity, but I calculate the value up such that the Integrand there becomes 10^-3 and is still falling from there, so I can basically put it to 0 from the point up and I can run my Integral from 0 (or here the value C_threshold = 10^-3 ) to up. The Parameter values used in my example are given by MassSquared = 100, Temp = 500 and Sign = +1. Is there Information missing which you need? — MrBaer, Jul 24 '15 at 08:45

Vivian Miranda · Accepted Answer · 2015-07-24T14:21:40.387

Whenever you ask a question that involves numerical results (and also bugs) you need to provide a concise code that we can analyze and run to reproduce your original problem. So yes, information is missing. We need

(1) concise code (we don't need the VEffClass!) with the GSL integration (you partially provided that - but your code is still quite complicated for stack overflow purposes).

(2) concise code with your implementation of the Simpson-rule.

(3) Snapshot of your maple calculation.

(2) and (3) are quite important because your question arises due to the fact that methods (2) and (3) disagree with (1), which you assumed is the incorrect answer. Why did you assume that (1) is incorrect and not the other way around? Do you have an analytical solution to check the numbers?

When I analyzed, in Wolfram Mathematica, the integrand with the parameters values you provided in a comment, I got an answer that is consistent with the GSL calculation. So I can't reproduce your problem and without (2) + (3) I can only guess what went wrong...

I attached a snapshot of my Mathematica notebook.

Thank you for the detailed reply! It was no Programm error, but through your reply and after trying to reproduce it, I found that there was a pre-factor of Temp^4/(2*Pi^2) missing between the two routines. — MrBaer, Jul 27 '15 at 07:49

GSL Incorrect values with numerical Integration

1 Answers1