9
$\begingroup$

I'm solving the Schrödinger equation for a particle in an annular geometry with hard wall boundary conditions and I've reduced it to the following equation:

$$J_m(k\,R_1)\,Y_m(k\,R_2) - J_m(k\,R_2)\,Y_m(k\,R_1) = 0,$$

where $R_1,R_2$ are the outer and inner radii respectively, and $J_m(r)$ and $Y_m(r)$ are the regular and irregular Bessel functions, labelled by integers $m$.

I want to find a way to create a function that is similar to BesselJZero, which returns the $n^{th}$ solution ($k_n$) of the above equation labelled by $m$. I tried using NSolve but I wasn't really able to construct something that works.

Just for future reference. I stumbled upon the following BesselZeros package. Apparently up until v.5.2 of Mathematica, they had a package called BesselZeros with functions such as BesselJYJYZeros[m,f,n] (where f = $R_{2}/R_{1}$) which returned the first $n$ zeros of precisely the function I was asking about.

As far as I can tell, this works well for $m\in \mathbb{Z}$ but runs into issues for complex values of the Bessel function index.

|improve this question
$\endgroup$
9
$\begingroup$

I've worked with a very similar equation quantization condition. Basically this is a root-finding problem, but since there are many roots and you want exactly the $n^\textrm{th}$, FindRoot is not very helpful. The approach which worked for me is to use NDSolve along with Sow and Reap. This question spells out the basics.

I'd try something like this:

quantcond[m_Integer, k_, R1_, R2_] := BesselJ[m, k R1] BesselY[m, k R2] - BesselJ[m, k R2] BesselY[m, k R1]

quantcondZero[m_?NumericQ, R1_?NumericQ, R2_?NumericQ, n_Integer, opts : OptionsPattern@NDSolve] := Module[
  {r = 0,
   min = 1/(1000 Min[R1, R2]),
   max = n 1000/Min[R1, R2]},

  Reap[
    NDSolve[
     {f'[k] == D[quantcond[m, k, R1, R2], k],
      f[min] == quantcond[m, min, R1, R2],
      WhenEvent[f[k] == 0, Sow[k]; r++],
      WhenEvent[r > n, "StopIntegration"]},
     f,
     {k, min, max},
     opts
    ]
   ][[2, 1]]
]

The function quantcondZero will return the first n zeros of the function quantcond.

Example

quantcondZero[1, 20^-3, 30^-3, 9]

{3160.94, 6293.06, 9431.39, 12571.3, 15711.9, 18852.9, 21994., 25135.2, 28276.5}

The roots look almost like multiples of π, which is the expected behavior.

Notes

  1. The main difficulty lies in scaling. Here, assuming that $R_1$ and $R_2$ are roughly the same scale, your zeros appear to scale roughly inversely with the size of $R_1$ and $R_2$, hence the setting for min and max, the limits of the integration. Depending on the behavior of the function, these may need to be adjusted to ensure you capture all of the zeros.
|improve this answer
$\endgroup$
7
$\begingroup$

I typically solve these problems by harnessing the power of Plot. Taking

fn[m_, k_, R1_, R2_] := BesselJ[m, k R1] BesselY[m, k R2] - BesselY[m, k R1] BesselJ[m, k R2]

we first Plot the function as

p1 = Plot[fn[1, k, 1, 2], {k, 0, 50}
  , PlotRange -> All
  , MeshFunctions -> {#2 &}
  , MeshStyle -> {Red, PointSize[0.02]}
  , Mesh -> {{0.}}]

yielding

enter image description here

We extract the roots and refine them using FindRoot:

pts = Sort@Cases[Normal@p1, Point[a_] :> First@a, Infinity]
k /. FindRoot[fn[1, k, 1, 2] == 0, {k, #}] & /@ pts

enter image description here

We can wrap this all up into a Module with the following code (borrowing Virgil's name from their answer):

quantcondZero[m_?NumericQ, R1_?NumericQ, R2_?NumericQ, n_Integer, x_?Positive] := Module[
  {plt, pts}
  , plt = Plot[BesselJ[m, k R1] BesselY[m, k R2] - BesselY[m, k R1] BesselJ[m, k R2]
      , {k, 0, x}, PlotRange -> All
      , MeshFunctions -> {#2 &}, Mesh -> {{0.}}]
  ; pts = Sort@Cases[Normal@plt, Point[a_] :> First@a, Infinity]
  ; k /. FindRoot[eqn[1, k, 1, 2] == 0, {k, #}] & /@ pts[[1 ;; Min[n, Length@pts]]]
  ]

Because one can't tell a priori how far out to plot the function, I have to include it as an argument to the function, and it makes it not completely automated. So, for instance, if we want the first 5 roots with $m=2$, we can start by doing

quantcondZero[2, 1, 2, 5, 10]
(* {3.40692, 6.42777, 9.52285} *)

It doesn't give us enough, so we try a larger x:

quantcondZero[2, 1, 2, 5, 20]
(* {3.40692, 6.42777, 9.52285, 12.6404, 15.7673} *)

Important note: this method is not the most robust. If you need a lot of zeroes, it is likely that the zeroes will get so close together that the Plot might miss them. You will need to do things like manually increase the number of PlotPoints in order to get it to work.

|improve this answer
$\endgroup$
6
$\begingroup$

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.

|improve this answer
$\endgroup$
  • $\begingroup$ For some reason, sometimes this function doubles the number of zeros. For example. annularZero[1, 2, 2, 2n+1] = annularZero[1, 2, 2, 2n] where $n\geq1$. But for $\rho_a = 0.1$ and $\rho_b = 2$, it won't do this anymore. I wasn't really able to understand what was causing this behaviour. $\endgroup$ – Aegon Jan 8 '16 at 1:00
  • $\begingroup$ @Aegon, I think I still have a factor wrong somewhere, so the bracketing still isn't quite right - will look at it some more... pretty sure it's fixable. $\endgroup$ – Jens Jan 8 '16 at 4:04
  • $\begingroup$ @Aegon Looks like the library add-on package also has trouble with duplicate roots. Maybe you can try my new approach and see if it works better now. $\endgroup$ – Jens Feb 1 '16 at 18:54
  • $\begingroup$ Thanks! This definitely works a lot better and so far seems to avoid duplicating and doesn't miss any roots. $\endgroup$ – Aegon Feb 3 '16 at 19:24

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.