[Top][All Lists]
[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
Re: [Bug-gsl] GSL interpolation
From: |
Jens Chluba |
Subject: |
Re: [Bug-gsl] GSL interpolation |
Date: |
Fri, 27 Jul 2012 04:12:13 -0400 |
On 2012-07-27, at 3:39 AM, Evgeny Kurbatov wrote:
> On Thu, 26 Jul 2012 14:39:21 -0400
> Jens Chluba <address@hidden> wrote:
>
>> Hi there,
>>
>> just a comment on the new convention for interpolation. Round-off
>> error lead to problems with
>>
>> if (x < interp->xmin || x > interp->xmax)
>> {
>> GSL_ERROR_VAL("interpolation error", GSL_EDOM, GSL_NAN);
>> }
>>
>> if indeed x very close to the boundaries but within the interval
>> [xmin, xmax]. It is probably better to put something like
>>
>> if (x < interp->xmin*(1.0-1.0e-14) || x >
>> interp->xmax*(1.0+1.0e-14)) {
>> GSL_ERROR_VAL("interpolation error", GSL_EDOM, GSL_NAN);
>> }
>>
>> This allows for a tiny range of extrapolation but fixes the round-off
>> error problem.
>
> Hi, Jens!
>
> I was the one who proposed a patch for interpolation routines
> toughening control to the boundary violation. There was times when I
> was sorry about that, but now the belief rises that this was right. The
> control forces you to think deeper and be more accurate.
> So hold on, be man!
It certainly is very important to make sure extrapolation is not mistaken for
interpolation. I absolutely agree that this is a good move! Still, [platform &
compiler flag (fast-math) dependent] round off errors should not be making you
have to dig into things. Certainly fixed this by manning up ;-)
>> In addition, I was wondering if the GSL team could be interested in
>> some pretty fast and accurate integration schemes based on the
>> Patterson quadrature rules? These are fully nested and converge very
>> rapidly. I have the required c++ code and will be happy to provide
>> them. Certainly, this can only be a base for GSL, since you guys have
>> much higher coding standards...
>
> I don't mind.
It is not my invention (just like Gauss-Kronrod is not yours), but it beats
almost any quadrature rule! NAG uses it, and it is one of the most efficient
algorithms I have come across. Just take a look. Even for me as an amateur
programmer it did not take a lot to beat your integration routines & those of
scipy in most cases [still getting precision to all digits]. It would be
foolish not to consider these impressive methods. I am sure you will find all
you need digging yourself then.
> Best witches,
> Evgeny Kurbatov
All the best & many thanks for your kind response,
Jens
-------------------------------------------------------------------------------------------------------------
Jens Chluba Tel. ++1 416 978 5946
or ++1 416 262 8224
Senior Research Associate Fax ++1 416 978 3921
CITA, University of Toronto .+'''+. .+'''+.
60 St. George Street mailto:address@hidden
Toronto, Ontario M5S 3H8 '+...+' '+...+'
Canada
http://www.cita.utoronto.ca/~jchluba/
-------------------------------------------------------------------------------------------------------------