:

## Brent's method for root finding

```Here is some standard alternative to Newton's method (and
thus may safe some homework ... so what).

It will find a root of f (I think it must be continuous C^1)
in the interval ax ... bx, if it has different signs at the
boundaries.

The code is more or less translated from netlib C library
(or similar).

Usage:

Digits:=16;

f:= x -> exp(x)-Pi;
zBrent( f, 0.0, 2.0);
'f(%)': '%'= evalf(%);

steps used = 10
1.144729885849400
f(1.144729885849400) = 0.

f:= x -> tan(x)-x;
eps:= 10.0^(-Digits+1);
zBrent( f, Pi/2 + eps, Pi/2 + Pi - eps);
'f(%)': '%' = evalf(%);

steps used = 11
4.493409457909064
-14
f(4.493409457909064) = -0.4 10

So it even works in ugly cases and starting only a little
bit off from boundary singularities.

##########################################################

zBrent:=proc(f, ax, bx)
local a,b,c, fa,fb,fc, cb, t1, t2, p,q,
prev_step, tol_act, new_step, iCounter, EPS, nSteps;

EPS:= 0;
nSteps:= 15*Digits;

a:= evalf(ax);
b:= evalf(bx);
fa:= evalf(f(a));
fb:= evalf(f(b));
c:= a;
fc:= fa;

if ( (fa < 0 and fb < 0) or (fa > 0 and fb > 0)) then
WARNING( cat(procname, ": root must be bracketed"));
return 'procname( args )'; #Exit Function
end if;

for iCounter from 1 to nSteps do
prev_step:= b - a;

if ( abs(fc) < abs(fb) ) then
a:= b;
b:= c;
c:= a;
fa:= fb;
fb:= fc;
fc:= fa;
end if;

tol_act := EPS * abs(b) + 0/3;
new_step:= (c - b) / 2;

if (abs(new_step) <= tol_act ) then
break
end if;

if (abs(prev_step) >= tol_act and abs(fa) > abs(fb)) then
cb:= c - b;
if (a = c) then
t1:= fb / fa;
p:= cb * t1;
q:= 1 - t1;
else
q:= fa / fc;
t1:= fb / fc;
t2:= fb / fa;
p:= t2 * (cb * q * (q - t1) - (b - a) * (t1 - 1));
q:= (q - 1) * (t1 - 1) * (t2 - 1);
end if;

if (p > 0) then
q:= -q;
else
p:= -p;
end if;

if (p < (0.75 * cb * q - abs(tol_act * q) / 2) and
p < abs(prev_step * q / 2)) then
new_step:= p / q;
end if;
end if;

if (abs(new_step) < tol_act) then
if (new_step > 0) then
new_step:= tol_act;
else
new_step:= -tol_act;
end if;
end if;

a:= b;
fa:= fb;
b:= b + new_step;
if( b = b + new_step ) then break; end if;
fb:= evalf(f(b));

if ((fb > 0 and fc > 0) or (fb < 0 and fc < 0)) then
c:= a;
fc:= fa;
end if

end do;

if nSteps <= iCounter then
WARNING(cat( 'procname', ": iteration limit has been reached"));
end if;
print(`steps used` = iCounter);

return b;
end proc;
``` ﻿