🚀 go-pugleaf

(*
This is a Mathematics V.2:0  programme allowing  the user to 
numerically compute all the (real) roots of any continuous 
function on a given interval. Faster versions, using Secant and
Newton's methods, are also available.  Complex case is in preparation.
Text copy of the paper can be obtained from: akyildiz@maths.oxford.ac.uk

Yilmaz Akyildiz

Beware:  RealInterval builtin function is not (yet) bug-free...
*)

*)
BeginPackage["IntervalBisectionMethod`"];

intbisec::usage =
"intbisec[f, left, right, tolerance] finds all the 
(real) roots of the function f on the closed 
interval [left, right] with accuracy up to the tolerance.  
It uses a version of interval bisection method. f must be 
provided in the form f = Function[variable, body]. 
For example,
intbisec[Function[x, x Sin[x] - 1/2], Pi, 20, 0.0001]
will find all the roots of  x * Sin[x] - 1/2 = 0 on the 
closed interval [Pi, 20] with accuracy correct up to 3 
decimal digits. The output will be displayed in terms of 
the Mathematica builtin function RealInterval[{a, b}] for 
the root which is contained in the closed interval [a, b]. 
We prefer the output in this format since Mathematica can 
do arithmetic with RealInterval objects.";

Begin["`Private`"];

lefthalf[RealInterval[{x_, y_}]]    := RealInterval[{x,(x+y)/2.}];             
righthalf[RealInterval[{x_, y_}]]   := RealInterval[{(x+y)/2.,y}];             
containszero[RealInterval[{x_,y_}]] :=  (x <= 0 <= y);                         
bisect[int_RealInterval]            := {lefthalf[int],righthalf[int]};         
                                                                               
                                                                               
intbisec[f_Function, a_, b_, tol_] :=                                         
   With[                                                                       
      {step =   Select[(bisect/@#)//Flatten, containszero[N[f[#]]]&]& },          
      Nest[step, {RealInterval[N[{a,b}]]}, Ceiling[Log[2,N[(b-a)]/tol]]]             
          ];                                                                   
End[];

EndPackage[]

(* end of  IntervalBisectionMethod *)