module numsolve; % root finding.
% Author: H. Melenk, ZIB, Berlin
% Copyright (c): ZIB Berlin 1991, all rights resrved
fluid '(!*noequiv accuracy!* !*exptexpand);
global '(iterations!* !*trnumeric erfg!*);
symbolic procedure rdsolveeval u;
% interface function;
begin scalar e,vars,y,z,p,r,erfg,mode,all,!*exptexpand;
integer n;
erfg:= erfg!*; erfg!* := nil;
u := for each x in u collect reval x;
u:=accuracycontrol(u,6,50);
if (all:='all memq u) then u:=delete('all,u);
e :=car u; u :=cdr u;
% construct the function vector.
e:=if eqcar(e,'list) then cdr e else list e;
e := for each f in e collect
if eqexpr (f:=reval f) then !*eqn2a f else f;
n := length e;
% starting point & variables.
if eqcar(car u,'list) then
u := for each x in cdar u collect reval x;
for each x in u do
<<if eqcar(x,'equal) then
<<z:=cadr x; y:=caddr x;
if eqcar(y,'!*interval!*) then mode:='i;
>> else <<z:=x; y:=random 100>>;
vars:=z . vars; p := y . p>>;
vars := reversip vars; p := reversip p;
if not(n=length vars) then
<<
% lprim "numbers of variables and functions don't match;";
% lprim "using steepest descent method instead of Newton";
% return rdsolvestdeval list ('list.e,'list.u,
% {'equal,'accuracy,accuracy!*},
% {'equal,'iterations,iterations!*});
rederr "numbers of variables and functions don't match"
>>;
if mode='i and length vars>1 or mode neq 'i and all then
rederr "mode for num_solve not implemented";
r:=if mode='i then rd!-solve!-interv(e,vars,p,n,all)
else rdnewton0(e,vars,p,n);
erfg!* := erfg;
return r;
end;
put('num_solve,'psopfn,'rdsolveeval);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% finding root in an interval (secant method and related)
%
symbolic procedure rd!-solve!-interv(e,vars,p,n,all);
(begin scalar u,fu,l,fl,x,acc,r,de,w,oldmode;
n := nil;
if length p>1 then typerr('list . p,"single interval");
p:=car p; e:=car e; x:=car vars;
w := {'boundseval,mkquote {e,{'equal,x,p}}};
if not memq(dmode!*,'(!:rd!: !:cr!:))then
<<oldmode:=t; setdmode('rounded,t)>>;
if errorp errorset(w,nil,nil) then
typerr(e,"bounded expression");
acc := !:!:quotient(1,expt(10,accuracy!*));
l:=evaluate(cadr p,nil,nil); u:=evaluate(caddr p,nil,nil);
fl:=evaluateuni(e,x,l); fu:=evaluateuni(e,x,u);
if not all then
dm!:(if zerop fl then return prepf l else
if zerop fu then return prepf u else
if fl*fu<0 then
r := regula!-falsi(e,x,l,fl,u,fu,acc));
if null r then de := reval {'df,e,x};
if null r and not all then
r:=rd!-solve!-trynewton(e,de,x,l,fl,u,fu,acc);
if null r then
r:=rd!-solve!-hardcase(e,de,x,l,fl,u,fu,acc,all);
if oldmode then setdmode('rounded,nil);
if r then return r;
rederr "not applicable";
end) where !*roundbf=!*roundbf;
symbolic procedure make!-rdsolve!-result1(x,u);
'list.for each r in u collect {'equal,x,prepf r};
symbolic procedure rd!-solve!-trynewton(e,de,x,l,fl,u,fu,acc);
begin scalar r,invde;
invde := {'quotient,1,de};
dm!: <<
if fu*evaluateuni(de,x,u) >0 then
r := rdnewton2({e},{{invde}},{x},acc,{u},'root,l,u);
if null r and fl*evaluateuni(de,x,l) <0 then
r := rdnewton2({e},{{invde}},{x},acc,{l},'root,l,u);
if r and (r:=car r) and l <= r and r <= u then
r := make!-rdsolve!-result1(x,{r}) else r:=nil;
>>;
return r;
end;
symbolic procedure regula!-falsi(e,x,l,fl,u,fu,acc);
make!-rdsolve!-result1(x,
{regula!-falsi1(e,x,l,fl,u,fu,acc,0,1)});
symbolic procedure regula!-falsi1(e,x,l,fl,u,fu,acc,mode1,mode2);
% modified regula falsi: using bisection in order to
% traverse root as often as possible.
dm!: begin scalar m,fm;
if mode1=mode2 then m:=(u+l)/2
else m := l*fu/(fu-fl) + u*fl/(fl-fu);
if (u-l) < acc then return m;
fm := evaluateuni(e,x,m);
if zerop fm then return m;
return if fl*fm<0
then regula!-falsi1(e,x,l,fl,m,fm,acc,nil,mode1)
else regula!-falsi1(e,x,m,fm,u,fu,acc,t,mode1);
end;
symbolic procedure mkminus u; {'minus,u};
symbolic procedure evaluateuni(e,x,p);
evaluate(e,{x},{p});
symbolic procedure dmboundsuni(e,x,l,u);
begin scalar i;
i:=boundseval {e,{'equal,x,{'!*interval!*,prepf l,prepf u} }};
return numr simp cadr i . numr simp caddr i;
end;
symbolic procedure rd!-solve!-hardcase(e,de,x,l,fl,u,fu,acc,all);
dm!: begin scalar il1,il2,pt,iv,r,rr,b,b1,q,d;
integer lev;
il1:={(l.fl) .(u.fu)};
main_loop:
lev:=lev+1; il2:=nil;
if null il1 then goto done;
il1 := reverse il1;
rd!-solve!-hardcase!-print(lev,il1);
loop:
if null il1 then goto bottom;
iv :=car il1; il1:= cdr il1;
l:=caar iv; fl:=cdar iv; u:=cadr iv; fu:=cddr iv;
if (d:=u-l) <0 then goto loop; %empty
if abs fl<acc then<<pt:=l;goto root_found>>;
if abs fu<acc then<<pt:=u;goto root_found>>;
b:=dmboundsuni(de,x,l,u);
% left boundary of interval
if (fl>0 and not((b1:=car b)<0))
or(fl<0 and not((b1:=cdr b)>0))
or not((pt:=l-fl/b1)<u) then goto loop;
if pt=l then goto precerr;
q:=evaluateuni(e,x,pt);
if not all and q*fl<0 then return regula!-falsi(e,x,l,fl,pt,q,acc);
if abs q<acc then goto root_found;
l:=pt; fl:=q;
% right boundary of interval
if (fu>0 and not((b1:=cdr b)>0))
or(fu<0 and not((b1:=car b)<0))
or not((pt:=u-fu/b1)>l) then goto loop;
if pt=u then goto precerr;
q:=evaluateuni(e,x,pt);
if not all and q*fu<0 then return regula!-falsi(e,x,pt,q,u,fu,acc);
if abs q<acc then goto root_found;
u:=pt; fu:=q;
% new point
pt :=(l+u)/2; %midpoint
q:=evaluateuni(e,x,pt);
if not all and q*fu<0 then return regula!-falsi(e,x,l,fl,pt,q,acc);
% generate new intervals
if not(abs q<acc) then
<<il2 :=find!-root!-addinterval(pt.q,u.fu,
find!-root!-addinterval(l.fl,pt.q,il2));
goto loop;
>>;
root_found:
r:=pt.r; if not all then goto done;
rr:=find!-root!-range(e,x,pt,acc);
il2 :=find!-root!-addinterval(cdr rr,u.fu ,
find!-root!-addinterval(l.fl,car rr,il2));
goto loop;
bottom:
il1 :=il2;
goto main_loop;
done:
r:=sort(r,function lessp);
return make!-rdsolve!-result1(x,r);
precerr:
rederr "precision not sufficient for finding all roots";
end;
symbolic procedure rd!-solve!-hardcase!-print(lev,il1);
if !*trnumeric then
<< prin2t {"recursion level:",lev,"remaining intervals:",length il1};
writepri( mkquote( 'list.
for each i in il1 collect
{'list,{'!*interval!*,prepf caar i,prepf cadr i},
dm!: prepf(cadr i-caar i)}),
'only);
>>;
symbolic procedure find!-root!-addinterval(l,h,il);
dm!: <<if car l < car h then il:=(l.h).il; il>>;
symbolic procedure find!-root!-range(e,x,p,acc);
% p is a point in (l .. u) where abs e(x)<acc.
% Find the next interval where e(x) > acc.
dm!:<<
while not(abs(fl:=evaluateuni(e,x,pl:=pl-acc/2))>acc) do;
while not(abs(fu:=evaluateuni(e,x,pu:=pu+acc/2))>acc) do;
(pl.fl).(pu.fu)
>> where pl=p,pu=p,fl=nil,fu:=nil;
if errorp errorset('(!:difference nil nil),nil,nil) then
<<
symbolic procedure !:difference(u,v);
if null u then !:minus v else if null v then u
else if u=v then nil
else if atom u and atom v then u-v
else dcombine(u,v,'difference);
symbolic procedure !:quotient(u,v);
if null u or u=0 then nil
else if null v or v=0 then rerror(poly,12,"Zero divisor")
else if atom u and atom v
% We might also check that remainder is zero in integer case.
then if null dmode!* then u/v else dcombine(u,!:recip v,'times)
else dcombine(u,v,'quotient);
>>;
endmodule;
end;