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

[Vipin K. Varma]

You're treating the wrong function in this last email unlike in the
codes you sent.

Vipin

On 10/05/2014 07:53 PM, Juan Pablo Amorocho wrote:
> Hi Vipin, 
> Taking a closer look these are my comments. If  in f(x) = log(cos(a) -
> cos(x)) / (x - a)  we have that a = pi/2, then f is not define between
> x \in[0, pi/2], and in general we have the condition cos(a) - cos(x) >
> 0 and x != a. So I really don’t think it is so much the method, but
> the the integration boundaries.  
> Here are a couple o links of QAGS and the quadrature it implements: 
>
>  http://www.netlib.org/quadpack/doc
>
> http://www.encyclopediaofmath.org/index.php/Gauss%E2%80%93Kronrod_quadrature_formula
> <http://www.encyclopediaofmath.org/index.php/Gauss%96Kronrod_quadrature_formula>
>
> Another point is the the error on the approximation depends on the
> high-order derivates of the function one is trying to integrate, and
> by judging the first derivate… oh boy! but Mathematica such help you
> as well.
>
> The point (if this were about that) goes to Klaus for pointing at the
> behaviour of the function. 
>
> I hope this helps, and if my argumentation is wrong, please (please!)
> let me know. 
>
> – Juan Pablo
> On 05 Oct 2014, at 17:53, Vipin K. Varma <address@hidden
> <mailto:address@hidden>> wrote:
>
>> 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.
>>>>>
>>>>>
>>>>
>>>
>>>
>>>
>>
>>
>