My goal was to give an answer using FindRoot
. In addition, I want to be able to uniquely identify the roots by the node number. I believe this is an important part of the original question, because it refers specifically to "labeling" the solutions. The method below achieves this unique association of the label n
with the number of zeros. For a given n
, it does not require finding all roots up to n
and counting them. This is where it differs fundamentally from the other answers.
The main issue is to find the starting points for FindRoot
reliably without duplicating or skipping roots. For large enough m
, this can be done very nicely by relying on BesselJZero
, but this doesn't work for all roots.
Therefore, I decided to modify my answer and not rely on BesselJZero
at all. The alternative is to use an asymptotic approach based on the WKB approximation. By doing this approximation with a rescaled coordinate $\xi \equiv \log(\rho)$, you get a very accurate formula for the starting points that works for large and small m
. I compared it to the solution available in the package linked in the question, and found it promising that my function doesn't duplicate roots whereas the library function does so:
BesselJYJYZeros[20, .5, 10]
(*
==> {25.4172, 25.4172, 29.9645, 34.0415, 38.1114, 42.5018, \
47.2888, 52.4004, 57.7525, 63.2819}
*)
You can see that the first root is duplicated, showing that the search was started with bad initial conditions.
Here is my approach:
Clear[annularZero, action, sWKB];
action[m_, k_, ρa_] :=
Sqrt[k^2 - m^2] - Sqrt[-m^2 + k^2 ρa^2] +
m ArcTan[(m (-Sqrt[k^2 - m^2] + Sqrt[-m^2 + k^2 ρa^2]))/(
m^2 + Sqrt[k^2 - m^2] Sqrt[-m^2 + k^2 ρa^2])]
sWKB[m_, k_?NumericQ, r_] :=
Re[action[m, k, Max[m/k, r]]];
Options[annularZero] = {"ShowPlot" -> False,
WorkingPrecision -> $MachinePrecision};
annularZero[m_, ρa_, n_, OptionsPattern[]] := Module[
{k, det, start1, start2, refine = False, firstPass, prec, showPlot},
(* For large m it helps to increase WorkingPrecision if FindRoot \
complains: *)
prec = OptionValue[WorkingPrecision];
(* With showPlot True, display the function whose roots we seek: *)
showPlot = OptionValue["ShowPlot"] =!= False;
(* Use WKB approximation to find two starting values for exact \
root search. No need for perfection,
so suppress warnings: *)
{start1, start2} =
Table[
Chop[k /.
Quiet@FindRoot[
sWKB[m, k, ρa] - Pi (n + δ), {k,
m/ρa}]], {δ, {-1/4, 0}}];
(* Want the n-th exact root of det: *)
det = SetPrecision[
BesselJ[m, k ρa] BesselY[m, k] -
BesselJ[m, k] BesselY[m, k ρa], prec];
(* Since det has infinitely many roots,
optionally show plot of interval around desired region: *)
If[showPlot,
Print[
Plot[det, {k, Max[0, start1 - 10], start2 + 10},
Epilog -> {Red, PointSize[.02],
Point[{{start1, 0}, {start2, 0}}]}, PlotRange -> {-1, 1}]]];
(* FindRoot with two starting values to bracket the desired root:*)
Quiet[
Check[
firstPass = Chop[k /. First@FindRoot[det, {k, start1, start2}]],
refine = True]];
(* If previous search issues warnings,
try to improve last result using Automatic method: *)
If [refine, If[showPlot, Print["refining zero"]];
Chop@N[k /.
First@FindRoot[det, {k, firstPass}, WorkingPrecision -> prec]],
firstPass]
];
For now, I simplified life for myself by requiring that the second argument ρa
must be less than 1. This can always be achieved by rescaling the original problem to appropriate length units. I may drop this restriction once the reliability has been tested some more. Here is the comparison to the test case above:
Table[annularZero[20, .5, i], {i, 10}]
(*
==> {25.4172, 29.9645, 34.0415, 38.1114, 42.5018, 47.2888, \
52.4004, 57.7525, 63.2819, 68.9443}
*)
Instead of duplicating the first root, I now get a new additional root at the end of the list. The other results are in complete agreement with the old add-on package.