3
$\begingroup$

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:

equation

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:

Missing zeroes

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:

all 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!

$\endgroup$
  • 1
    $\begingroup$ Did you try to use shifting border limits - it potentially can allow to resolve the lost roots? I mean something like series of scans with slightly different star/stop points and further join the results excluding duplicates.. $\endgroup$ – Rom38 Mar 16 '15 at 14:05
  • $\begingroup$ That did it! I had tried making the intervals shorter but in line with what you suggested I did a series of scans of length 0.01 for each $d_B$. I now find all the zeroes - the gaps are filled up. I will post the graphic above. It takes a long time to produce, but it's accurate. Thank you $\endgroup$ – Latrace Mar 16 '15 at 22:04
  • $\begingroup$ Scanning methods are the typical solution in optics when you have met the physical limit of a resolution.. $\endgroup$ – Rom38 Mar 17 '15 at 22:05
2
$\begingroup$

You were using an older version of the function I posted in the linked answer. With the new version, copied from that answer, it's possible to specify the option PlotPoints to increase the resolution of the plot grid that's used for bracketing. This in principle allows you really find all roots. The tradeoff is that the speed decreases with increasing PlotPoints. Since even the original code in the question runs very long, I only ramped the points up to 500 to show how it improves the plot in principle. Caution - this takes a few minutes to run:

For[Subscript[d, B] = end, Subscript[d, B] >= start, 
 Subscript[d, B] -= prec, 
 roots = Sort[findAllRoots[Detf[e, n], {e, -1, 0}, PlotPoints -> 500],
    Less];
 Do[PrependTo[LIST[[i]], {N[Subscript[d, B]], roots[[i]]}], {i, 
   Length[roots]}];]

LIST = DeleteCases[LIST, {}];
ListPlot[LIST, 
 AxesLabel -> {"\!\(\*SubscriptBox[\(d\), \(B\)]\)", "E"}]

roots

As you can see, the number of missing roots is much smaller with this additional option. Now it's up to you to increase it to 1000 if you have the patience.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.