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

[Vipin K. Varma]

Yes, it's indeed curious that x=0 seems to misbehave. One can avoid the
point x=0 altogether for the case a=pi/2 because then the function is
symmetric about x=pi/2; and so we can integrate from [pi/2, pi] and then
simply double the result i.e. perform the following:

  gsl_integration_qags (&F, M_PI/2, M_PI, 0, 1e-12, 1000, w, &result,
&error);
  printf ("result          = % .18f\n", 2*result);

 This gives -0.454682346139364646, which is correct to 14 digits
(compared with Mathematica); but why the need to avoid the x=0 point
(which is clearly equivalent to x=pi for a=pi/2) is unclear to me. I
cannot, however, afford to treat a=pi/2 as a special case; therefore any
explanations for the failure of QAGS is appreciated.

Vipin

On 10/05/2014 03:27 PM, Juan Pablo Amorocho D. wrote:
> Hi all,
>
> I'm baffled! If I plot the function I get an asymptote at x = 0, but 
> evaluating the function with the code I sent hours earlier I get the
> pair (x,f(x)) (0.000000, -0.456196)...
>
>  In any case, if the  integration range is 0 <x <= pi   and I replace
> 0 for 1e-6 as in
>
>   gsl_integration_qags (&F, 0, M_PI, 0, 1e-6, 1000, w, &result, &error);
>                                         ^ changed to 1e-6
>
> I get the value result          = -0.454681894556977995
> here is what I actually ran:
>
> #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;
>   int n = (params->n);
>   double a = (params->a);
>   double f = log((cos(a) - cos(x))/(x - a));
>   return f;
> }
>
> int main (void)
> {
>   gsl_integration_workspace * w
>     = gsl_integration_workspace_alloc (1000);
>   int n = 1;
>   size_t neval;
>   double result, result_int, error;
>   double alpha = M_PI/2;
>   struct my_f_params params = {n, alpha};
>
>   gsl_function F;
>   F.function = &f2;
>   F.params = &params;
>
>   gsl_integration_qags (&F, 1e-6, M_PI, 0, 1e-6, 1000, w, &result,
> &error);
>
>   printf ("result          = % .18f\n", result);
>
>   gsl_integration_workspace_free (w);
>
>   return 0;
> }
>
>
>
> Comments are highly appreciate!
>
> -- Juan Pablo
>
> 2014-10-05 13:32 GMT+02:00 Vipin K. Varma <address@hidden
> <mailto:address@hidden>>:
>
>     Hi all,
>
>     Thanks very much for the replies.
>
>     @Klaus: The argument to the logarithm is zero only when either
>     a=-x or a=x=0,pi; neither situation is possible because 0<a<pi,
>     while 0<=x<=pi for the integration range.
>
>     @Juan Pablo: Indeed I have tried larger relative errors but the
>     integration always fails when a = 1*M_PI/2; and as your code
>     shows, there is no undefined function value close to, or equal to,
>     that value of 'a'. This is my code without any indentation [the
>     cos(n*x) part of the function can be ignored]:
>
>     <CODE>
>     #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;
>       int n = (params->n);
>       double a = (params->a);
>       double f = /*cos(n*x)**/log((cos(a) - cos(x))/(x - a));
>       return f;
>     }
>
>     int main (void)
>     {
>       gsl_integration_workspace * w
>         = gsl_integration_workspace_alloc (1000);
>
>       int n(1);
>       size_t neval;
>       double result, result_int, error;
>       double alpha = M_PI/2;
>       struct my_f_params params = {n, alpha};
>
>       gsl_function F;
>       F.function = &f2;
>       F.params = &params;
>
>       gsl_integration_qags (&F, 0, M_PI, 0, 1e-6, 1000, w, &result,
>     &error);
>
>       printf ("result          = % .18f\n", result);
>
>       gsl_integration_workspace_free (w);
>
>       return 0;
>     }
>     <CODE>
>
>     Best,
>     Vipin
>
>
>     On 10/05/2014 08:53 AM, Juan Pablo Amorocho wrote:
>>     HI all, 
>>     I evaluated the function on the integration, and didn’t see
>>     anything suspicious. Though one thing caught my eye. Vipin, you
>>     are passing  a relative error of 10e-12. Have you tried a value
>>     of 10e-6 or 10e-7? By the way, what’s is the value of n and what
>>     does int n(1) resolves to? In the code below I assume n = 1.  I
>>     don’t mean to be picky, but could you please send us (or me)
>>      code that we could “just” copy& paste without further
>>     modification in order to compile and run? Thanks!
>>
>>     — Juan Pablo
>>
>>     #include<stdio.h>
>>     #include<math.h>
>>
>>
>>     doublef (doublex){
>>     const double a = 0.992*M_PI/2;
>>     constdoublen = 1.0;
>>     return cos(n*x)*log((cos(a) - cos(x))/(x - a));
>>
>>     }
>>
>>     int main(){
>>
>>     double a = 0.0;
>>     double h = 0.01;
>>     double fa = 0.;
>>     do{
>>     fa = f(a);
>>     printf("(%f, %f)\n", a, fa);
>>     a+=h;
>>     }while(a <=M_PI);
>>
>>     return0;
>>
>>     }
>>
>>     On 04 Oct 2014, at 19:08, Klaus Huthmacher
>>     <address@hidden
>>     <mailto:address@hidden>> wrote:
>>
>>>     Dear Vipin,
>>>
>>>>     //  double f = cos(n*x)*log((cos(a) - cos(x))/(x - a));//
>>>
>>>     Is it possible that the argument of your logarithm can be zero
>>>     for PI? Or
>>>     that you divide by zero, if x-a becomes zero for PI?
>>>
>>>     Kind regards,
>>>      Klaus.
>>>
>>>
>>
>
>
>