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 = ¶ms; > > 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 <vva...@ictp.it > <mailto:vva...@ictp.it>>: > > 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 = ¶ms; > > 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 >> <huthmac...@physik.uni-kl.de >> <mailto:huthmac...@physik.uni-kl.de>> 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. >>> >>> >> > > >