AllRoots
- To: MATHGROUP at yoda.physics.unc.edu
- Subject: AllRoots
- From: "Mathematical Institute, (0865) 2-73525" <AKYILDIZ at vax.ox.ac.uk>
- Date: Thu, 8 Oct 92 15:15 BST
(*
Thnaks to Troels Petersen, here is a slightly better version.
There has been a good number of interest in finding ALL the
roots of any equation by the interval bisection method.
Below is a program which does the job, (sorry for the
indentation, I am never a good programmer perhaps someone
can help me on that). The text of the paper can still be
obtained via e-mail from akyildiz at vax.oxford.ac.uk
*)
(*
This program allows the user to numerically compute
all the (real) roots of any continuous function on
a given interval.
*)
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], 0, 20, 0.0001]
will find all the roots of x * Sin[x] - 1/2 = 0 on the
closed interval [0, 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`"];
rght[RealInterval[{x_, y_}]] := y;
middle[RealInterval[{x_, y_}]] := (x + y)/2.;
lft[RealInterval[{x_, y_}]] := x;
width[RealInterval[{x_, y_}]] := y - x;
includes[RealInterval[{a_, b_}], RealInterval[{c_, d_}]] :=
a <= c && d <= b;
contains[RealInterval[{a_, b_}], c_] := a <= c && c <= b;
intbisec[f_Function, left_, right_, tolerance_] :=
Module[{e = RealInterval[{-tolerance, tolerance}],
j, k = 1, xL, xR, x, y, zeroes = {}},
x[1] = RealInterval[{left, right}];
While[ k != 0,j = 0;
Do[If[width[x[i]] < tolerance ||
includes[e, x[i]],
zeroes = Append[zeroes, x[i]],
xL[i] = RealInterval[{lft[x[i]], middle[x[i]]}];
xR[i] = RealInterval[{middle[x[i]], rght[x[i]]}];
If[contains[f[xL[i]], 0], j++; y[j] = xL[i]];
If[contains[f[xR[i]], 0], j++; y[j] = xR[i]]
], {i, k}];
k = j; Do[x[n] = y[n], {n, k}]
]; zeroes
]
End[];
EndPackage[]
(* end of IntervalBisectionMethod *)