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

[Juan Pablo Amorocho D.]

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.
>
>
>
>
>
> --
> ==============================================================
> 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: address@hidden
> ==============================================================
>
>