Re: [Gretl-devel] c code question

On 16 Oct 2017, at 17:35, Allin Cottrell <cottrell(a)wfu.edu&gt; wrote: On Mon, 16 Oct 2017, Riccardo (Jack) Lucchetti wrote: > [T]he quality of the result is comparable to what we already get via Richardson extrapolation. The complex-step technique performs a second-order approximation, which is what we get by setting fdjac_quality to 2. For example, try the following script (which uses the definition of rderiv as per Oleh's previous script): > > <hansl> > set verbose off > > function scalar cube(matrix *x) > return x[1]^3 > end function > > g = rderiv("x^3", 0) > printf "rderiv: %36.33f\n", g > > x = {0} > > set fdjac_quality 1 > g = fdjac(x, cube(&x)) > printf "fdjac(1): %36.33f\n", g > > set fdjac_quality 2 > g = fdjac(x, cube(&x)) > printf "fdjac(2): %36.33f\n", g > </hansl> > > On my system I get > <output> > rderiv: -0.000000000000000000000000000000049 > fdjac(1): 0.000000000000000888178419700125430 > fdjac(2): 0.000000000000000000000000000000000 > </output> OK, but try Oleh's rderiv() examples: <hansl> d1 = rderiv("sin(x)",$pi/6) - cos($pi/6) x = {$pi/6} d2 = fdjac(x,(sin(x))) - cos($pi/6) printf "sin(x): d1 = %g, d2 = %g\n", d1, d2 d1 = rderiv("exp(2*x)",10) - 2*exp(20) x = {10} d2 = fdjac(x,exp(2*x)) - 2*exp(20) printf "exp(2x): d1 = %g, d2 = %g\n", d1, d2 d1 = rderiv("x+x^2+x^3",1) - (1 + 2 + 3) x = {1} d2 = fdjac(x,(x+x^2+x^3)) - (1 + 2 + 3) printf "cubic: d1 = %g, d2 = %g\n", d1, d2 d1 = rderiv("exp(x)",100) - exp(100) x = {100} f = fdjac(x, exp(x)) d2 = f - exp(100) printf "exp(100): d1 = %g, d2 = %g\n", d1, d2 printf " (proportional error = %g)\n", abs(f - exp(100)) / exp(100) d1 = rderiv("100+x",0.1) - 1 x = {0.1} d2 = fdjac(x,(100+x)) - 1 printf "100+x: d1 = %g, d2 = %g\n", d1, d2 d1 = rderiv("sin(x)",(20000+1/4)*$pi) - cos((20000+1/4)*$pi) x = {(20000+1/4)*$pi} d2 = fdjac(x, sin(x)) - cos((20000+1/4)*$pi) printf "sin(x): d1 = %g, d2 = %g\n", d1, d2 </hansl> With fdjac_quality 2 I get: <output> sin(x): d1 = 0, d2 = -2.88153e-11 exp(2x): d1 = 0, d2 = -0.0899086 cubic: d1 = 0, d2 = -2.70636e-10 exp(100): d1 = 0, d2 = -4.09506e+32 (proportional error = 1.52339e-11) 100+x: d1 = 0, d2 = -3.15016e-09 sin(x): d1 = 0, d2 = -1.72044e-07 </output> There's no doubt that the complex method gives greater accuracy. However, as you say, this method is limited to functions that accept complex arguments and produce complex return values.

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