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

[Vipin K. Varma]

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