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

[Juan Pablo Amorocho]

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

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> 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>:
>> 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>
>>> 
>>> 
>>> double f (double x){
>>> const double a = 0.992*M_PI/2;
>>> const double n = 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);
>>> 
>>> return 0;
>>> 
>>> }
>>> 
>>> On 04 Oct 2014, at 19:08, Klaus Huthmacher <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.
>>>> 
>>>> 
>>> 
>> 
>> 
>> 
> 
>