I am interested in bracketing the first zero
(starting somewhere) of an arbitrarily
complicated function as a preliminary to
using find_root or bf_find_root.
The traditional method is shown here,
using the change in sign of the function,
but I need another method (see below)
in addition.
I adapt Numerical Recipes code for zbrak,
stopping with the first root found:
syntax: zbrak (expr,x,x1,x2,num)
returns list [xr1,xr2] of positions which
bracket the first root found, working from
x1 to x2, or else the empty list [].
num is the number of intervals the
domain [x1, x2] is divided into.
In floating point arithmetic, basing a decision
on the sign of the product of a function evaluated
at two different, but close, values may be
risky if the two points are close to the
zero of the function, leading to underflow of
the product. Hence better to just compare
signs of f(x1) and f(x2) rather than looking
at the sign of f(x1)*f(x2); hence the
function opp_sign. (Trying to use Maxima's
sign function is awkward.)
=================================
Maxima 5.28.0-2 http://maxima.sourceforge.net
using Lisp GNU Common Lisp (GCL) GCL 2.6.8 (a.k.a. GCL)
(%i1) opp_sign (o1,o2) :=
(if o1 < 0 and o2 > 0 then true
else if o1 > 0 and o2 < 0 then true
else false)$
(%i2) zbrak(ze,zx,zx1,zx2,zn1) :=
block([zdx,zfp,zxp,zfc,zr:[],zn ],local(zf),
define(zf(zx),ze),
if evenp(zn1) then zn : zn1 + 1
else zn : zn1,
zxp : zx1,
zdx : (zx2-zx1)/zn,
zfp : zf(zxp),
for zi:1 thru zn do
(zxp : zxp + zdx,
zfc : zf(zxp),
if opp_sign (zfp,zfc) then
(zr : [zxp - zdx,zxp],
return()),
zfp : zfc),
zr)$
(%i3) zbrak(x,x,-1,1,3);
(%o3) [-1/3,1/3]
(%i4) zbrak(sin(x),x,-1,1,21);
(%o4) [-1/21,1/21]
(%i5) zbrak(sin(x),x,-1,1,20);
(%o5) [-1/21,1/21]
================================
This approach will not work at bracketing
a root of a function which has, for example,
f(x1) = 0, f'(x1) = 0, and f''(x1) > 0
(for example, concave upwards at x = 1,
zero slope at x = 1, and having a local
minimum at x = 1) such as f(x) = (x - 1)^2 .
Another example which will not work with the
above approach is f(x) = 10/cos(10*x)^2
I need an approach to the latter two cases
which doesn't require picking apart the function,
but just basing the method of bracketing the root on
general behaviors of the function, thus allowing
the function to be arbitrarily complicated
(including the presence of special functions).
The motivation is improving the nint package
(Maxima by Example, Ch. 8) by looking for
possible non-integrable singularities of the
integrand inside the quadrature domain.
(Checking for the presence of endpoint sigularities
and their integrability is easier and already
worked out.)
My initial approach for interior singularities
is to look for zeros of the inverse of the
first derivative of the given integrand.
Any suggestions would be appreciated.
Ted Woollett