Dear Johan,
I think you are correct. The intention was to use an absolute error that is just
a dummy so that the relative error is what is enforced by the integration
routine. The code was written under the (as it turns out incorrect) assumption
that the gsl_integration_qag routine would guarantee the tighter of the two
given limits, but instead it may return if either one of them is fulfilled. So
one should indeed pass 0 to make the absolute error limit a dummy...
Fortunately, in this particular case the integrated function is so smooth that
the default 41 Gauss-Kronrod already returned a sufficiently accurate result, so
that this issue did not cause a problem for a gadget2. But it should
nevertheless be changed of course (which I in fact did for gadget3).
Thanks for pointing this out,
Volker
Johan Maes wrote:
> Hm apparently the approximation is ok if either one of the error limits
> is satisfied:
> http://linux.math.tifr.res.in/manuals/html/gsl-ref-html/gsl-ref_16.html#SEC244
>
> /Each algorithm computes an approximation to a definite integral of the
> form, /
>
> /I = \int_a^b f(x) w(x) dx
> /
>
> / where w(x) is a weight function (for general integrands w(x)=1). The
> user provides absolute and relative error bounds (epsabs, epsrel) which
> specify the following accuracy requirement, /
>
> /|RESULT - I| <= max(epsabs, epsrel |I|)
> /
>
> / where RESULT is the numerical approximation obtained by the algorithm.
>
> /So, in order to make the absolute error a dummy, one could just set it
> to 0....right?
>
> Johan
> /
> /Johan Maes wrote:
>> Hi all,
>>
>> This may be a minor thing, well this is a minor thing, but still, I'm
>> wondering why in init_drift_table (in driftfac.c), the Hubble constant
>> is used as absolute error for the integrations. Here's a description
>> of adaptive gsl integration taken from
>> http://linux.math.tifr.res.in/manuals/html/gsl-ref-html/gsl-ref_16.html#SEC246
>>
>> /This function applies an integration rule adaptively until an
>> estimate of the integral of f over (a,b) is achieved within the
>> desired absolute and relative error limits, epsabs and epsrel
>> /
>> In the code it says the absolute error is just used as a dummy, so I
>> guess the goal is to make it as big as possible (compared to the
>> result) so the above is always ok. But since the results have
>> dimension of time, I would rather expect some fraction of the Hubble
>> time then...what am I missing here? I'm doing something similar now to
>> convert a to t, using the integration in Gadget as an example, that's
>> why I noticed.
>>
>> Thx & cheers,
>>
>> Johan
>> --
>> En toen zei de kikker: "Voor mij ne kleine me stoverij, alstublieft."
>>
>
> --
> En toen zei de kikker: "Voor mij ne kleine me stoverij, alstublieft."
>
Received on 2009-03-20 22:26:37