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

[Vipin K. Varma]

Hi Juan Pablo,

Thanks for the inputs, indeed better precision may be obtained by
treating the x=a condition separately in the function definition.

Best,
Vipin

On 10/08/2014 04:13 AM, Juan Pablo Amorocho D. wrote:
> 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)
>


--