Re: [Bug-gsl] GSL interpolation

[Jens Chluba]

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/
-------------------------------------------------------------------------------------------------------------