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'm making a wrapper

for quadtable to do easy integration

with hansl, and my simplistic sfd2*

works systematically better, while

it should not to.

Oleh

*) fd from fdjac

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

> Dear all, First, a couple of numbers to introduce the question
> cos(0.1)~ 0.99500416 5278 025766095561987803870294838576225415084035959352...
> ### NOTE all results are on 32-bit Ubuntu 15.10 To compare a helper:
> function string as_string (scalar x) sprintf s "%.15g",x return s end function
> 1)  x={0.1} eval as_string(fdjac(x,sin(x))) 0.99500416 7780 28
> 2)
> set fdjac_quality 1 eval as_string(fdjac(x,sin(x))) 0.99500416 5451 974
> 3)
> set fdjac_quality 2 eval as_string(fdjac(x,sin(x))) 0.99500416 3123 667
> A home-backed thing inspired by fdjac help :
> A helper: function scalar feval (string fun, scalar x) return @fun end function
> The main:
> function scalar sdj2 (string f, scalar x) if abs(x)>10^3 h = 10^3*sqrt(x/1000*$macheps) else h=10^(3)*sqrt($macheps) endif z = (8*(feval(f,x+h)-feval(f,x-h)) - (feval(f,x+2*h)-feval(f,x-2*h)))/(12*h) sz = as_string(z) return z end function
> We have
> eval as_string(sdj2("sin(x)",0.1)) 0.99500416 5277 274 recall: 0.99500416 5278 025766
>
> In theory, fdjac should be better, since it uses the same formula as sdj2() and it should optimize the step h, but... Oleh

In general, numerical differentiation should never be expected to be 
accurate to more than a few decimal digits, for the combined effect of the 
inherent imprecision of machine precision and discrete (as opposed to 
_truly_ infinitesimal) differentiation. Run this script and you'll see 
what I mean:

<hansl>
set echo off
set messages off
clear

set fdjac_quality 0
loop for (x=-1; x<1; x+=2/13)
     matrix a = {x}
     printf "x = %5.2f, precision = %d\n", a, \
     floor(-log10(abs(1-fdjac(a, sin(a))/cos(x))))
endloop

set fdjac_quality 1
loop for (x=-1; x<1; x+=2/13)
     matrix a = {x}
     printf "x = %5.2f, precision = %d\n", a, \
     floor(-log10(abs(1-fdjac(a, sin(a))/cos(x))))
endloop

set fdjac_quality 2
loop for (x=-1; x<1; x+=2/13)
     matrix a = {x}
     printf "x = %5.2f, precision = %d\n", a, \
     floor(-log10(abs(1-fdjac(a, sin(a))/cos(x))))
endloop
</hansl>

That said, I would have expected Richardson extrapolation to fare somewhat 
better in this case. I'll investigate.

-------------------------------------------------------
   Riccardo (Jack) Lucchetti
   Dipartimento di Scienze Economiche e Sociali (DiSES)

   Università Politecnica delle Marche
   (formerly known as Università di Ancona)

   r.lucchetti@univpm.it
   http://www2.econ.univpm.it/servizi/hpp/lucchetti
-------------------------------------------------------

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