Writing a bisection root finder

Question

I'm trying to make a simple function that will calculate roots of equations by using bisection method. My code isn't working right. Can anyone take a look and help me?

Here is the code:

findroot[funk_Symbol, {var_ , init_, a_, b_, xi_}, eps_] := 
  While[Abs[funk[init]] >= eps;
    xi = N @ (a + b)/2;
    If[funk[xi] > 0, a = a; b = xi, a = xi, b = b];
    {init -> xi};
    Abs[funk[xi]]]
f[x_] := x^2 - 3;
findroot[f, {x, 3, 1, 2, xi}, 0.00001]

varibles are:

funk_Symbol -> function defined by user,
var_ -> variable (x in this case),
init_ -> initial guess,
a_ and b_ -> left and right border,
xi -> (a+b)/2
eps -> error defined by user

Related: Implement the Bisection algorithm elegantly and easily — C. E., Apr 09 '18 at 19:52

m_goldberg · Answer 1 · 2018-04-09T20:38:10.827

This is both simpler and corrects some errors in your formulation.

Clear[findroot]
findroot[
    func : (_Symbol | _Function),
    {a_?NumericQ, b_?NumericQ},
    eps_?Positive] :=
  Module[{x, y, aa = N @ a, bb = N @ b},
    While[Abs[y = func[x = (aa + bb)/2]] >= eps,
      If[Sign[y] === Sign[func[aa]], aa = x, bb = x]];
    x]

Passing a named function as the 1st argument

f[x_] := x^2 - 3
findroot[f, {1, 2}, 1.*^-8]

1.73205

Passing a pure function as the 1st argument

findroot[Cos[#] - # &, {.5, 1}, 1.*^-8]

0.739085

Notes

One reason why your code doesn't work is that in Mathematica you can not assign to the the formal arguments of a function in its code body because those names don't exist in the code body when it is evaluated. The values of the actual arguments passed are substituted for the formal arguments before the code body evaluation begins.
Another reason is that you must test for x being on same side of the root as aa. That condition holds when the sign of func[x] is the same as the sign of func[aa].

@AccidentalFourierTransform. Good point. Now that I think about it, I don't like it myself. I have made the change you suggest. — m_goldberg, Apr 09 '18 at 20:39
Just to add: it's more robust to use x = aa + (bb - aa)/2 (per this one). — J. M.'s missing motivation, Apr 10 '18 at 02:13

score 2 · Accepted Answer · edited Apr 10 '18 at 04:32

Since it seems to be an academic exercise, here's an alternative in a different style using patterns and replacements:

ClearAll[findroot, init, left, right, done];

init[f_] := {
   {a_, 0, b_, _} :> a,   (* done *)
   {a_, _, b_, 0} :> b,   (* done *)
   {a_, sa_, b_, sb_} :> ({a, sa, b, sb, #, Sign@f[#]} &[(a + b)/2])};
left[f_] := {a_, sa_, b_, sb_, c_, sc_} /; sa*sc == -1 :>
   ({a, sa, c, sc, #, Sign@f[#]} &[(a + c)/2]);
right[f_] := {a_, sa_, b_, sb_, c_, sc_} /; sb*sc == -1 :>  
   ({c, sc, b, sb, #, Sign@f[#]} &[(b + c)/2]);
done[] = {a_, _, b_, _, c_, sc_} /; a === b || sc == 0 :> c;

findroot[func : (_Symbol | _Function), {a_, b_}, eps_: $MachineEpsilon] /; eps > 0 := 
  With[{obj = Chop[func[#], eps] &},
   {a, Sign@obj@a, b, Sign@obj@b} /. init[obj] //. {done[], left[obj], right[obj]}
   ];

Note: It has some checks missing in the OP. For instance, init[] checks whether the initial values a or b are roots; done[] stops when the interval {a, b} cannot be subdivided further (well, up to Internal`$SameQTolerance). Furthermore, if func at the initial {a, sa, b, sb, c, sc} has all the same signs at a/b/c, it will return the data structure showing the same signs. (One could catch this and throw an error.) Otherwise, after done[] checks whether c is a solution, the replacement rule left[] subdivides the left half, and right[] subdivides the right half. I used Chop via obj to implement the tolerance eps. It also never re-evaluates the function at the same point.

Examples:

findroot[Cos, {1., 2.}]  (* uses default eps = $MachineEpsilon *)
% - Pi/2
(*
  1.5707963267948966` 
  0.
*)

findroot[Cos, {1., 2.}, 0.001]
% - Pi/2
(*
  1.5703125` 
  -0.000483827
*)

f[x_] := x^2 - 3;
findroot[f, {1., 3.}, 0.00001]
% - Sqrt[3]
(*
  1.73205
  -2.7729*10^-6
*)

Writing a bisection root finder

2 Answers2