Re: [Gretl-devel] fdjac: a better example

On Tue, 29 Dec 2015, oleg_komashko(a)ukr.net wrote: Yes, with the status-quo small step size in fdjac2.c, plus Richardson, the results can become violently unstable for small values of the function argument. Further investigation: <hansl> set echo off set messages off function scalar fbig (scalar y) return y + 10000 end function scalar f0 f1 f2 loop for (x=1.0; x>1e-10; x/=10) -q y = {x} set fdjac_quality 0 f0 = fdjac(y,fbig(y)) set fdjac_quality 1 f1 = fdjac(y,fbig(y)) set fdjac_quality 2 f2 = fdjac(y,fbig(y)) printf "y=%.2e: %f, %f, %f\n", x, f0, f1, f2 endloop </hansl> Here, this gives <gretl-output> y=1.00e+00: 1.000000, 1.000000, 1.000000 y=1.00e-01: 0.999756, 1.000061, 0.999858 y=1.00e-02: 1.000977, 1.000977, 1.001994 y=1.00e-03: 0.976563, 0.976563, 0.996908 y=1.00e-04: 1.220703, 0.915527, 0.101725 y=1.00e-05: 0.000000, 0.000000, 2.034505 y=1.00e-06: 0.000000, 0.000000, 40.690104 y=1.00e-07: 0.000000, 0.000000, 0.000000 y=1.00e-08: 0.000000, 0.000000, -2034.505208 y=1.00e-09: 0.000000, 0.000000, 20345.052083 y=1.00e-10: 0.000000, 0.000000, 101725.260417 </gretl-output> It's unsurprising that the "basic" fdjac methods give an answer of zero when y gets small enough: y + h + 10000 becomes numerically indistinguishable from y + 10000 when you only have 64 bits per floating-point number. But the Richardson output goes crazy in a way that's not very nice at all. In current git I've modified fdjac2 to use a step size as in your script (when Richardson is selected). We need to do more testing on this in case it goes nasty under some conditions, but here's what I'm now getting from the script above: <gretl-output> y=1.00e+00: 1.000000, 1.000000, 1.000000 y=1.00e-01: 0.999756, 1.000061, 1.000000 y=1.00e-02: 1.000977, 1.000977, 1.000000 y=1.00e-03: 0.976563, 0.976563, 1.000000 y=1.00e-04: 1.220703, 0.915527, 1.000000 y=1.00e-05: 0.000000, 0.000000, 1.000000 y=1.00e-06: 0.000000, 0.000000, 1.000000 y=1.00e-07: 0.000000, 0.000000, 1.000000 y=1.00e-08: 0.000000, 0.000000, 1.000000 y=1.00e-09: 0.000000, 0.000000, 1.000000 y=1.00e-10: 0.000000, 0.000000, 1.000000 </gretl-output> Allin