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 <vva...@ictp.it > <mailto:vva...@ictp.it>> 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 = ¶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. >>>>> >>>>> >>>> >>> >>> >>> >> >> >
