Finding all the roots of a function of one real variable has been discussed extensively in e.g. Find all roots of an interpolating function (solution to a differential equation) (answer no. 2 by Jens provides a function findAllRoots that works without modification for some simple functions I tried). I'm dealing with quite a complicated function, however, and I am certain that findAllRoots is not finding all the roots. I have also tried the well-known RootSearch.m package and Stan Wagon's FindAllCrossings[] (e.g. About multi-root search in Mathematica for transcendental equations answer no. 3 by J. M.), but to no avail.
My general question is if anyone else has had a missing roots problem before, and, if so, how you solved it. I've tried playing with the WorkingPrecision AccuracyGoal PrecisionGoal and MaxIterations options in the FindRoot function but again this made no difference.
To be concrete I should mention the function who's roots I'm after:
Here $T_N$ and $U_{N-1}$ are the Chebyshev polynomials, $d_B$ is a parameter to be varied in the interval $[0,5]$ and $E$ is the unknown. I am interested in this function because it appears in a published physics article, this one, but background in physics is not necessary to solve the problem I'm having, I believe. Unless that the region of interest is $E \in (-\infty,0)$, but with the parameter choice I mention $E \in (-1,0)$. Here is my code:
Clear[findAllRoots]
SyntaxInformation[
findAllRoots] = {"LocalVariables" -> {"Plot", {2, 2}},
"ArgumentsPattern" -> {_, _, OptionsPattern[]}};
SetAttributes[findAllRoots, HoldAll];
Options[findAllRoots] =
Join[{"ShowPlot" -> False, PlotRange -> All},
FilterRules[Options[Plot], Except[PlotRange]]];
findAllRoots[fn_, {l_, lmin_, lmax_}, opts : OptionsPattern[]] :=
Module[{pl, p, x, localFunction, brackets},
localFunction = ReleaseHold[Hold[fn] /. l :> x];
If[lmin != lmax,
pl = Plot[localFunction, {x, lmin, lmax},
Evaluate@
FilterRules[Join[{opts}, Options[findAllRoots]],
Options[Plot]]];
p = Cases[pl, Line[{x__}] :> x, Infinity];
If[OptionValue["ShowPlot"],
Print[Show[pl, PlotLabel -> "Finding roots for this function",
ImageSize -> Large, BaseStyle -> {FontSize -> 8}]]], p = {}];
brackets =
Map[First,
Select[(*This Split trick pretends that two points on the curve \
are "equal" if the function values have _opposite _ sign.Pairs of \
such sign-changes form the brackets for the subsequent FindRoot*)
Split[p, Sign[Last[#2]] == -Sign[Last[#1]] &],
Length[#1] == 2 &], {2}];
x /. Apply[FindRoot[localFunction == 0, {x, ##1}] &,
brackets, {1}] /. x -> {}];
Subscript[k, A][e_] = Sqrt[Subscript[V, 0] + e];
Subscript[k, B][e_] = Sqrt[-e];
Detf[e_, n_] :=
Exp[-n*Subscript[k, B][e]*Subscript[d,
B]] (ChebyshevT[n,
Cos[Subscript[k, A][e] Subscript[d, A]] Cosh[
Subscript[k, B][e] Subscript[d, B]] - (
Subscript[k, A][e]^2 - Subscript[k, B][e]^2)/(
2 Subscript[k, A][e] Subscript[k, B][e])
Sin[Subscript[k, A][e] Subscript[d, A]] Sinh[
Subscript[k, B][e] Subscript[d, B]]] - (Cos[
Subscript[k, A][e] Subscript[d, A]] Sinh[
Subscript[k, B][e] Subscript[d, B]] - (
Subscript[k, A][e]^2 - Subscript[k, B][e]^2)/(
2 Subscript[k, A][e] Subscript[k, B][e])
Sin[Subscript[k, A][e] Subscript[d, A]] Cosh[
Subscript[k, B][e] Subscript[d, B]])*
ChebyshevU[n - 1,
Cos[Subscript[k, A][e] Subscript[d, A]] Cosh[
Subscript[k, B][e] Subscript[d, B]] - (
Subscript[k, A][e]^2 - Subscript[k, B][e]^2)/(
2 Subscript[k, A][e] Subscript[k, B][e])
Sin[Subscript[k, A][e] Subscript[d, A]] Sinh[
Subscript[k, B][e] Subscript[d, B]]]);
prec = 10^-2;
start = 0;
end = 5;
Subscript[d, A] = 1;
Subscript[V, 0] = 1;
n = 30;
LIST = Table[{}, {50}];
For[Subscript[d, B] = end, Subscript[d, B] >= start,
Subscript[d, B] -= prec,
roots = Sort[findAllRoots[Detf[e, n], {e, -1, 0}], Less];
Do[PrependTo[LIST[[i]], {N[Subscript[d, B]], roots[[i]]}], {i,
Length[roots]}];
]
LIST = DeleteCases[LIST, {}];
ListPlot[LIST,
AxesLabel -> {"\!\(\*SubscriptBox[\(d\), \(B\)]\)", "E"}]
(sorry about the subscript notation, it looks better in the notebook). This code searches for the roots $E$ for a given $d_B \in [0,5]$ in increments of 0.01 using findAllRoots. It's easy to use another algorithm, just replace findAllRoots[...] in the roots = Sort[findAllRoots[Detf[e, n], {e, -1, 0}], Less]; line by your root-searching function of choice. I get the following error message when running this:
FindRoot::brmp: The root has been bracketed as closely as possible with machine precision but the function value exceeds the absolute tolerance 1.0536712127723497`*^-8.
but I don't know what to do about it. The figure the code produces is this:
You see the gaping blank holes of missing zeroes! (yes - they are missing). Any help would be greatly appreciated! Thank you.
A solution was provided by Rom38: for each $d_B$, do a set of scans on small $E$ intervals and throw away duplicates in the end. This gave me the sought-after zeroes:
The bulk of the code is the same, and I replaced the bottleneck by this:
For[Subscript[d, B] = end, Subscript[d, B] >= start,
Subscript[d, B] -= prec,
For[j = -1.1, j < 0, j += prec,
roots =
Sort[findAllRoots[Detf[e, n], {e, j, j + prec}], Less];
Do[PrependTo[LIST[[i]], {N[Subscript[d, B]], roots[[i]]}], {i,
Length[roots]}];
]
]
LIST = DeleteCases[LIST, {}];
LIST = DeleteDuplicates[LIST];
I aborted the calculation before $d_B$ got to zero but it's more than enough to see that it works. Thank you Rom38!
