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>: > 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> > > > 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 <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. > > > > > > -- > ============================================================== > Vipin Kerala Varma > > Postdoctoral fellow, > Condensed Matter and Statistical Physics Section, > The Abdus Salam International Centre for Theoretical Physics, > Strada Costiera 11, > Trieste, Italy 34151 > > Phone: +39-040-2240-369 > > email: vva...@ictp.it > ============================================================== > >