[Gretl-devel] Another extra.gfn candidate?
Sorry for the bombardment: I just realised that, much to my surprise, not
only we don't have a numerical root-finding function in hansl; it seems we
don't have one in libgretl either!
So I put together a quick hansl prototype for the crudest algorithm there
is (the bisection method). It's reportd below, with a few examples. The
question is very much similar to the one in the other thread: shall we put
this into extra.gfn, or add this (hopefully, in a refined form) to
libgretl and eventually expose a user-level function?
<hansl>
set verbose off
function scalar bisection_solve(string f, scalar a[0], scalar b[NA])
if ok(b)
# initialisation via a bracket
top = xmax(a,b)
bot = xmin(a,b)
x = (top+bot)/2
else
# initialisation via a guess
x = a
w = abs(x) .* 1.1 + 1
bot = x - w
top = x + w
endif
scalar iter = 0
scalar MAXITER = 100
scalar conv = 0
scalar eps = 1.0e-13
loop while !conv && iter<MAXITER --quiet
x0 = x
scalar y = feval(f, x)
printf "Iter %4d: [%8.3f,%8.3f] f(%g) = %g\n",
iter++, bot, top, x, y
if y * feval(f, top) < 0
bot = x
else
top = x
endif
x = (top+bot)/2
conv = abs(x - x0) < eps
endloop
return x
end function
### example with default initial guess
function scalar f(scalar x)
return 3 * x^3 - 1
end function
x = bisection_solve("f")
printf "\nf(%14.11f) = %g (exact = %14.11f)\n\n\n", x, f(x), 3^(-1/3)
### example with given initial guess
function scalar f(scalar x)
return log(x) + x
end function
x = bisection_solve("f", 0.5)
printf "\nf(%14.11f) = %g (exact = %14.11f)\n\n\n", x, f(x), \
0.56714329040978387299996866221 # so-called "omega" constant
### example with given bracket
function scalar f(scalar x)
return sin(x) - 1/abs(1+cos(x))
end function
x = bisection_solve("f", -1, 1)
printf "\nf(%14.11f) = %g\n\n\n", x, f(x)
</hansl>
-------------------------------------------------------
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
-------------------------------------------------------