[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
Re: [Help-gsl] Numerical integration QAGS cannot reach tolerance and abo
From: |
Vipin K. Varma |
Subject: |
Re: [Help-gsl] Numerical integration QAGS cannot reach tolerance and aborts |
Date: |
Sun, 05 Oct 2014 17:53:08 +0200 |
User-agent: |
Mozilla/5.0 (X11; Linux x86_64; rv:31.0) Gecko/20100101 Thunderbird/31.1.2 |
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 <address@hidden
> <mailto: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 = ¶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
>> <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.
>>>
>>>
>>
>
>
>
- [Help-gsl] Numerical integration QAGS cannot reach tolerance and aborts, Vipin K. Varma, 2014/10/04
- Re: [Help-gsl] Numerical integration QAGS cannot reach tolerance and aborts, Klaus Huthmacher, 2014/10/04
- Re: [Help-gsl] Numerical integration QAGS cannot reach tolerance and aborts, Juan Pablo Amorocho, 2014/10/05
- Re: [Help-gsl] Numerical integration QAGS cannot reach tolerance and aborts, Vipin K. Varma, 2014/10/05
- Re: [Help-gsl] Numerical integration QAGS cannot reach tolerance and aborts, Juan Pablo Amorocho D., 2014/10/05
- Re: [Help-gsl] Numerical integration QAGS cannot reach tolerance and aborts,
Vipin K. Varma <=
- Re: [Help-gsl] Numerical integration QAGS cannot reach tolerance and aborts, Juan Pablo Amorocho, 2014/10/05
- Re: [Help-gsl] Numerical integration QAGS cannot reach tolerance and aborts, Vipin K. Varma, 2014/10/06
- Re: [Help-gsl] Numerical integration QAGS cannot reach tolerance and aborts, Juan Pablo Amorocho D., 2014/10/07
- Re: [Help-gsl] Numerical integration QAGS cannot reach tolerance and aborts, Vipin K. Varma, 2014/10/08