# This procedure finds a zero of a function using bisection. The parameters
# are an expression in an unknown defining the function, the name of the
# unknown, a range known to contain a zero, and a tolerance value for the
# accuracy to which the root is found. The zero is found using floating
# point computation.
#
# The value of the function at one of the end-points of the range passed
# must be less than or equal to zero, and the other must be greater than
# or equal to zero.
#
# If the expression defining the function contains a call of a procedure,
# it will probably be necessary to put it in single-quotes (') to prevent
# immediate evaluation.
#
# If a very small tolerance value is passed as the last parameter, the
# zero will be computed to full floating-point accuracy, except if it
# very close to zero.
#
# Example: bisection (x^2-2, x, 0..10, 1e-100) should return sqrt(2).
bisection := proc(f,x,rng,tolerance)
local lx, hx, mx, dx, lf, hf, mf;
if tolerance<=0 then
ERROR(`Fourth operand must be a positive tolerance value`);
fi;
if not type(rng,range) then
ERROR(`Third operand must be a range`);
fi;
lx := evalf(op(1,rng));
hx := evalf(op(2,rng));
if hx0 and hf>0 or lf<0 and hf<0 then
ERROR(`Range does not straddle a zero of the function`);
fi;
do # loop until we find zero, or interval is smaller than tolerance
dx := hx-lx;
mx := lx + dx/2;
if abs(mx-lx)0 and lf>0 or mf<0 and lf<0 then
lx := mx;
else
hx := mx;
fi;
od;
end: