Re: [Gretl-devel] c code question
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