For defining h  I used the following

if abs(x)>10^3
        h = 10^3*abs(x)/1000*sqrt($macheps)

else
h=10^(3)*sqrt($macheps)

1000 is more or less recommended

by the literature.

With |x|>1000 it will be the same

I thought, dfjac uses several

calculations to define h.

Interesting. what, say, lapack does use?

Oleh

On Tue, 29 Dec 2015, oleg_komashko@ukr.net wrote:

> Dear Riccardo, My simple script uses exactly Richardson 
> extrapolation but it uses a simplistic method for defining the 
> step. 3 digits better is by far too large. Clearly, fdjac uses a 
> more advanced way to find h. It looks, as if it has a buglet 
> inside.

I don't think there's any bug in Jack's Richardson code. In fact I 
believe the only reason you're getting more accurate results in this 
case is that you got lucky with your "simplistic method" for setting 
h. If I just change the code in our minpack-derived fdjac2.c to use 
your choice of

h = 1000 * sqrt(DBL_EPSILON)

in place of

h = fabs(x[j]) * sqrt(DBL_EPSILON)

I get a value for sin'(0.1) which is good to 12 decimal places!

I don't suppose this would be consistent across arbitrary functions, 
but it may be worth experimenting with a bigger h, or a different 
formula for h, in the "quality = 2" case.

Allin
_______________________________________________
Gretl-devel mailing list
Gretl-devel@lists.wfu.edu
http://lists.wfu.edu/mailman/listinfo/gretl-devel