Re: [Gretl-devel] fdjac precision question
On Tue, 29 Dec 2015, oleg_komashko(a)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(a)univpm.it
http://www2.econ.univpm.it/servizi/hpp/lucchetti
-------------------------------------------------------