[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 14:16:08 -0400

On 2012-07-27, at 4:56 AM, Evgeny Kurbatov wrote:

>>> 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 ;-)
> You're right on the one hand but the tolerance value you're propose
> fixates precision on some poorly controlled level. Why 1e-14, not
> GSL_DBL_EPSILON? On the other hand, this tolerance may cause necessity
> to make the range check in routines where the interpolation results are
> used. Namely such a hardly to catch bug has led to the patch.
> It's just my opinion.

Maybe I did misunderstand, but previously for x < xmin the function returned 
f(xmin) and for x>xmax it gave f(xmax). That of course can lead to some 
difficulties, but the version I propose is just getting rid of possible cases 
where x==xmin and x==xmax is just not fulfilled because of round-off errors 
(which depend on compliers, optimization flags etc, if I understand correctly). 
If you prefer it of course can be  

  if (x < interp->xmin*(1.0-GSL_DBL_EPSILON) || x > 
      GSL_ERROR_VAL("interpolation error", GSL_EDOM, GSL_NAN);

But I am sure you guys even have a better solution for this issue, which I 
believe should be addressed because many users probably are surprised if they 
use an interpolation routine and then get an error when they call it exactly at 
the boundary (which is not a very unnatural thing to do)!

>>>> 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.
> As for me I'm not a maintainer of GSL and didn't contribute to the
> integration routines.
> It would be easier for all, if you provide a code preferentially in a
> plain C because of the architecture of GSL. But also send a C++
> version, please.

Below my c++ implementation. As I said, I am an amateur programmer, but the 
code works very well. The main feature is that the quadrature rules a fully 
nested and not even one function evaluation is wasted until the maximum 
refinement level is reached. The quadrature rules are explained in Patterson T. 
N. L., 1968, Math. Comput., 22, 847.

>>> Best witches,
>>> Evgeny Kurbatov
>> All the best & many thanks for your kind response,
>> Jens

All the best,


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                        '+...+'          '+...+'


reply via email to

[Prev in Thread] Current Thread [Next in Thread]