Re: [Help-gsl] Numerical integration QAGS cannot reach tolerance and aborts

[Juan Pablo Amorocho D.]

Vipin,

sorry for the late reply and pointing out my mistake. I fiddle a bit more
and my conclusion is that there is nothing wrong this the method itself. It
just can't deliver the requested error conditions. This code and results
explains my point. The answer to why the quadrature performs so poorly lies
on the high-order derivative of the function.

I hope this helps!

-- Juan Pablo
#include <stdio.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_integration.h>

struct my_f_params { int n; double a; };

double f2 (double x, void * p) {
  struct my_f_params * params = (struct my_f_params *)p;
  double a = (params->a);
  double f = log((cos(a) - cos(x))/(x - a));
  return f;
}

int main (void)
{
  size_t neval;
  double result, result_int, error;
  double alpha = M_PI_2;
  struct my_f_params params = {alpha};
  gsl_integration_workspace * w
    = gsl_integration_workspace_alloc (1000);

  gsl_function F;
  F.function = &f2;
  F.params = &params;
  double errs[5] = {1e-2, 1e-3, 1e-4, 1e-6, 1e-7};
  double abserr[5] = {1e-1, 1e-2, 1e-3, 1e-5, 1e-8};

for(int i = 0; i < 5; i++){
  gsl_integration_qag (&F, 0, M_PI, errs[i], abserr[i], 1000, 6, w,
&result, &error);

  printf ("err = %e \tabserrs = %e \tresult = % .18f\t error= %f \n",
errs[i], abserr[i], result, error);
}
  gsl_integration_workspace_free (w);

  return 0;
}


gcc -std=c99 -L/usr/local/lib code_a.c -lgsl -lm -lgslcblas && ./a.out

err = 1.000000e-02     abserrs = 1.000000e-01     result =
-2.631950213795384741     error= 0.164315
err = 1.000000e-03     abserrs = 1.000000e-02     result =
-2.632228658613667172     error= 0.018008
err = 1.000000e-04     abserrs = 1.000000e-03     result =
-2.632263464212754034     error= 0.002247
err = 1.000000e-06     abserrs = 1.000000e-05     result =
-2.632268397653268366     error= 0.000018
gsl: qag.c:261: ERROR: could not integrate function
Default GSL error handler invoked.
Aborted (core dumped)