Bisection

Andy Long, Spring, 2024

Includes a little  debugging code, and verbose printing of the output if desired.

This is the printout command for verbose results:

Printout[output_,a_,b_,c_,fc_,slope_]:=
Module[{},
Print[NumberForm[TableForm[output,
TableHeadings->{None,{"(Subscript[a, k])","(Subscript[c, k])","(Subscript[b, k])","f[(Subscript[c, k])]","slope"}}],16]]; 
Print["  c  = ",NumberForm[c,16] ]; 
Print[" Δc  = ±",((b-a)/2)]; 
Print["f[c] = ",NumberForm[fc,16] ]; 
Print["f'[c] = ",NumberForm[slope,16] ]; 
]

Bisection code:

(* These are options that one can give at game time: 
-- epsilon: the value we consider "zero" (0 by default)
  -- tolerance: the width of the interval we consider "close enough" to the root
		-- verbose: print the details of the convergence
-- debug: print out a little more for the coder....
*)
Options[Bisection]={maxits->100,epsilon->0,tolerance->10^(-6),verbose->True,debug->False} ;

(* Here's the function itself: *)
Bisection[f_,a0_,b0_,OptionsPattern[]] :=

Module[{a,b,c,fa,fb,fc,i=1,its=OptionValue[maxits],eps=OptionValue[epsilon],tol=OptionValue[tolerance],verbose=OptionValue[verbose],debug=OptionValue[debug]},

(* use local variables (so as not to alter passed arguments *)
If[b0<a0,
{a,b}={b0,a0},
{a,b}={a0,b0}
];
{fa,fb} = {f[a],f[b]};  
If[debug,
Print[DecimalForm[{{a,fa},{b,fb}},20]]
];

(* Check first to see if there's a root at an endpoint: *)
If[Abs[fa]<=eps, 
If[debug,
Print[DecimalForm[fa,20]];
Print[DecimalForm[eps,20]];
];
If[verbose,
Print["Root at endpoint a."];
Return[Transpose[{{a,a,b,0}}]],
Return[a]
]
];
If[Abs[fb]≤eps,
If[debug,
Print[DecimalForm[a,20]];
Print[DecimalForm[fa,20]];
Print[DecimalForm[b,20]];
Print[DecimalForm[fb,20]];
Print["ArcSin[b]"];
Print[b>=Pi/2];
Print[DecimalForm[ArcSin[b],20]];
Print[DecimalForm[eps,20]];
];
If[verbose,
Print["Root at endpoint b."];
Return[Transpose[{{a,b,b,0}}]],
Return[b]
]
];

(* Check that the endpoints produce function values of opposite sign: *)
If[fa*fb>=0,Print["Error: function value doesn't change sign."];Return[Transpose[{{a0,a0,b0,0}}]]];

(* Here's our initial interval width: *)
width=b-a;
slope=(fb-fa)/width;

(* We don't need to do any MORE iterations than required by the tolerance, or by the user: *)
its = Min[its, Ceiling[Log2[width/tol]]];

(* c is the midpoint: *)
c = a+0.5*width; 
fc=f[c]; 

(* Let's collect the output, since we've gone to the trouble of computing it: *)
output={{a,c,b,fc,slope}}; 

(* Did we stomp on a root by good fortune? *)
If[Abs[fc]<=eps, 
If[verbose,
Print["Root at c."];Printout[output,a,b,c,fc,slope];Return[Transpose[{{a,c,b,fc}}]],
Return[c]
]
];

(* Otherwise, let's start looping: *)
While[a < c && c<b  && i<=its && width>tol,
If[ Sign[fa] ==  Sign[fc], 
a = c;fa=fc, 
b = c;fb=fc
];
width=width/2;
slope =(fb-fa)/width;
i=i+1;

(* Compute new midpoint: *)
c = a+0.5*width;
fc=f[c];

(* Stomp on root? *)
If[Abs[fc]≤eps, 
If[verbose,
Print["Root at c."];
Printout[output,a,b,c,fc,slope];
Return[Transpose[{{a,c,b,fc}}]],
Return[c]
]
];
(* Otherwise, store output and do it again: *)
output=Append[output,{a,c,b,fc,slope}]; 
]; 

(* We've popped out, so print some nice formatted output: *)
If[verbose,
Printout[output,a,b,c,fc,slope];
(* Return the output, so that it can be plotted, etc. *)
Return[Transpose[output]],
Return[c]
];
]

We first illustrate that the method works (at least on this problem!):

f[x_]:=Cos[x]-x
results=Bisection[f,0,Pi/2,maxits->10];
ListPlot[{results[[2]],results[[4]]},PlotStyle->{Thick,PointSize[Large]}]

What if our inputs don’t have function values of opposite sign?

results=Bisection[f,Pi/3,Pi/4,maxits->5];
ListPlot[{results[[2]],results[[4]]},PlotStyle->{Thick,PointSize[Large]}]

Now show that it can be used to find quotients, using only multiplication, addition, and subtraction:

f[x_]:=3x-1
results=Bisection[f,.1,1];
ListPlot[{results[[2]],results[[4]]},PlotStyle->{Thick,PointSize[Large]}]

Our interval endpoints are roots, and now we catch them:

f[x_]:=Sin[x]
results=Bisection[f,0,Pi];
ListPlot[{results[[2]],results[[4]]},PlotStyle->{Thick,PointSize[Large]}]
results=Bisection[f,0.1,Pi];
ListPlot[{results[[2]],results[[4]]},PlotStyle->{Thick,PointSize[Large]}]

Our first step is onto a root, but we miss it -- because of machine error! Our algorithm converts the answers to decimal, which has the effect of ruining a perfectly good root. Notice how long it takes us to get even somewhat close to the initial value of f[c] (our final result isn’t as good as our first calculated approximation).

f[x_]:=Cos[x]
results=Bisection[f,0,Pi];
ListPlot[{results[[2]],results[[4]]},PlotStyle->{Thick,PointSize[Large]}]

We can catch it by reducing our demand on f[c]: we give an epsilon of the appropriate size that we desire:

f[x_]:=Cos[x]
results=Bisection[f,0,Pi,epsilon->10^-6];
ListPlot[{results[[2]],results[[4]]},PlotStyle->{Thick,PointSize[Large]}]

What if we try to find a root of a discontinuous function?

f[x_]:=.5-Floor[x]
results=Bisection[f,-.1,Pi,epsilon->10^-6];
ListPlot[{results[[2]],results[[4]]},PlotStyle->{Thick,PointSize[Large]}]

Preparing to create a fake arsine, using our sinner fake of sine, and the interval [0,Pi/2]:

f[x_]:=Sin[x]-.5
x0=0;
x1=Pi/2;
results=Bisection[f,x0,x1,epsilon->10^-6,verbose->False];
results

x0=.999999
f[x_]:=Sin[x]-x0
Bisection[f,x0,x1,epsilon->10^-9,verbose->True];
ArcSin[x0]

x0=.99999
f[x_]:=Sin[x]-x0
Bisection[f,x0,x1,epsilon->10^-9,tolerance->10^-6,maxits->1000,verbose->False,debug->True];
ArcSin[x0]

What if we screw up, and write the endpoints in the wrong order? Bisection switches them...

Bisection[f,x1,x0,epsilon->10^-9,tolerance->10^-8,maxits->10,verbose->False,debug->True];