- Mail actions: [ respond to this message ] [ mail a new topic ]
- Contemporary messages sorted: [ by date ] [ by thread ] [ by subject ] [ by author ] [ by messages with attachments ]

From: Volker Springel <volker_at_MPA-Garching.MPG.DE>

Date: Fri, 20 Mar 2009 22:26:36 +0100

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

Date: Fri, 20 Mar 2009 22:26:36 +0100

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:

Received on 2009-03-20 22:26:37

*
This archive was generated by hypermail 2.3.0
: 2023-01-10 10:01:30 CET
*