I'm trying to find the roots of a function involving Bessel functions.

Here is my code

Clear[f, eta, chi, omega, alpha, beta, h, M, T, v, rho, phi, a, n, k, Zeros, z];
eta = Sqrt[beta^2 (Sqrt[4 alpha^2 omega^2/beta^4 + 1] - 1)/(2 alpha^2)];
chi = Sqrt[beta^2 (Sqrt[4 alpha^2 omega^2/beta^4 + 1] + 1)/(2 alpha^2)];
beta = Sqrt[alpha^2*(14.682 phi)/a^2];
alpha = a^2*omega/k;
a = .5;
f[n_, k_, phi_] = eta*(BesselJ[n - 1, eta*a] - BesselJ[n + 1, eta*a]) - 
   chi*((BesselI[n - 1, chi*a] + BesselI[n + 1, chi*a])*(BesselJ[n,eta*a]))/BesselI[n, chi*a];

Plotting the functions for different values of $n$ and $\phi$

Manipulate[Plot[f[n, k, phi], {k, 1, 300}, PlotRange -> {-10^1.5, 10^1.5}], 
     {n,0, 5, 1}, {phi, .25, 3, .25}]

enter image description here

Now, based on this post FindRoot giving false roots with Bessel Functions , I tried to obtain a list of the roots just by renaming my function and manually introduce the values of $n$ and $\phi$ as follows

n = 0;
phi = .25;
k = z;
Zeros[z_] = f[n, k, phi];
(zeros = First@
  Quiet@NDSolve[{y'[x] == Zeros'[x], y[0] == Zeros[0], 
     WhenEvent[y[x] == 0, Sow[FindRoot[Zeros[z], {z, x}]]]}, 
    y, {x, 0, 40}, AccuracyGoal -> 1, 
    PrecisionGoal -> 1]) // AbsoluteTiming

But it is not working!!! Mathematica gives the following error

First::first: {} has a length of zero and no first element. >>
{0.044039, First[{}]}

Any support that you can give me is highly appreciated.

  • 1
    $\begingroup$ Exactly what are the values of n and phi when this happens? And anything else that will allow someone to get to exactly where you are finding this problem. That will allow someone to try to reproduce your FindRoot not finding any roots, thus your Reap being empty, thus your First of empty list failing. When something doesn't work then stripping off the Manipulate and all the things you have added to try to hide the details of the failures may help diagnose exactly what the problem is. $\endgroup$ – Bill Sep 6 '14 at 15:02

I just used this as another test case for a different approach that I use a lot: it's based on my answer to Find all roots of an interpolating function.

First copy the definition of findAllRoots from that answer, then just do this:


roots {11.3873,41.2847,90.7297,159.864,248.719}

The option "ShowPlot"->True can be omitted if you don't want to se the plot displayed.

This shows the complete list of roots in the given interval {x, 0, 300}. Regarding the original source of the code for findAllRoots, I noted what I recalled in the linked answer, but I never was able to locate the old reference.

  • $\begingroup$ Sweet solution! +1 both of your sols/ans. Thanks a bunch! $\endgroup$ – Your Majesty Jun 24 '15 at 15:39

For example with n = 2.5 and phi = 0.5:

sol = FindInstance[f[2.5, k, .5] == 0 && 10 < k < 300, k, Reals, 15] // Values // Flatten

{45.6831, 100.726, 175.08, 269.028}

p = Point[Transpose[{sol, Table[0, {Length @ sol}]}]];

Plot[f[2.5, k, .5], {k, 1, 300}, Epilog -> {Red, PointSize[0.02], p}]

enter image description here

  • 1
    $\begingroup$ (+1) because the question doesn't ask specifically for all the roots. But if you do want all the roots in an interval, FindInstance doesn't always work. For example, with a small modification to f[0, k, .25], this same approach yields multiple approximate copies of the same roots, but misses the one at 248.719. $\endgroup$ – Jens Sep 6 '14 at 16:46
  • $\begingroup$ Thanks eldo! But this method is rather slow compared to the one proposed by Jens. $\endgroup$ – Jemme Sep 8 '14 at 22:09

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.