Posted on 2008-08-31 15:35 By
Axel Vogt (
987)
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;
1 hour 35 min ago
1 hour 55 min ago
3 hours 6 min ago
6 hours 11 min ago
6 hours 18 min ago
11 hours 49 min ago
19 hours 8 min ago
19 hours 32 min ago
19 hours 42 min ago
21 hours 44 min ago