Solve: $j_n(ka)=\frac{ka}{n}j_{n-1}(ka)$. Zeros of equation with spherical Bessel functions.

I am trying to solve for the boundary conditions of a resonant spherical cavity (with the centre cut out of it) and got to the point where I have boundary conditions that I am trying to solve:

$$$$j_n(ka)=\frac{ka}{n}j_{n-1}(ka),$$$$ and $$$$j_n(kb)=\frac{kb}{n}j_{n-1}(kb),$$$$

where $$j_n(x)$$ is the spherical Bessel function of the first kind, $$k$$ is the wave number, $$a$$ is the inner radius of the cavity, $$b$$ is the outer radius of the cavity, and $$n$$ (an integer) is the order of the Bessel function. Note that the original condition was $$j_n(ka)=-\frac{ka}{n}j_{n}'(ka)$$, where the derivative is w.r.t. $$a$$.

The spherical Bessel function is defined as:

$$j_n(x)=\sqrt{\frac{\pi}{2x}}J_{n+1/2}(x),$$

where $$J_n(x)$$ is the Bessel function of the first kind.

I am basically trying to solve the above two boundary conditions for $$k$$ given $$a$$ and $$b$$ (or perhaps $$a/b$$) or solve for $$a$$ and $$b$$ given $$k$$ (for fixed n). A solution to each condition individually (solve for $$ka$$ and $$kb$$) would work too.

Attempt #1

When I divide one equation by the other I get:

$$\frac{a}{b}=\frac{j_n(ka)}{j_n(kb)}\frac{j_{n-1}(kb)}{j_{n-1}(kb)},$$

but I don't know how to isolate for any of the variables.

Attempt #2

I was also looking at other definitions for spherical Bessel functions:

$$j_n(x)=(-x)^{n}\left(\frac{d}{x~dx}\right)^n\frac{\sin(x)}{x},$$

by substituting this definition into the boundary condition I get (note the funny $$x=ka$$ substitution):

$$(-ka)^{n}\left(\frac{d}{k^2a~da}\right)^n\frac{\sin(ka)}{ka}=\frac{ka}{n}(-ka)^{n-1}\left(\frac{d}{k^2a~da}\right)^{n-1}\frac{\sin(ka)}{ka}$$

Canceling the $$(-ka)$$ terms and factoring out a derivative operator:

$$-\left(\frac{d}{k^2a~da}\right)\left(\frac{d}{k^2a~da}\right)^{n-1}\frac{\sin(ka)}{ka}=\frac{1}{n}\left(\frac{d}{k^2a~da}\right)^{n-1}\frac{\sin(ka)}{ka}.$$

Now if I define $$f(a)=\left(\frac{d}{k^2a~da}\right)^{n-1}\frac{\sin(ka)}{ka}$$, then:

$$\frac{df}{da}=-\frac{k^2a}{n}f,$$

which is a separable equation (?) that can be solved by integrating:

$$\int\frac{1}{f}df=-\frac{k^2}{n}\int a da,$$

to get:

$$f=e^{-\frac{k^2a^2}{2n}},$$

which is interesting that I get a Gaussian, but I am not sure it is helpful (or even correct).

Attempt #3

I was also looking at recursion formulae as well as orthogonality conditions, but just ended up going in circles. Here is what I found with the first recursion formula:

$$j_n(ka)=\frac{ka}{2n+1}(j_{n-1}(ka)+j_{n+1}(ka)$$

substitute in the boundary condition for $$j_{n-1}(ka)$$:

$$j_n(ka)=\frac{ka}{2n+1}(\frac{n}{ka}j_n(ka)+j_{n+1}(ka))$$

now if I shift the index in the boundary condition up by one (this part suspicious):

$$j_{n+1}=\frac{ka}{n+1}j_n(ka),$$

which came be put into the above equation to give:

$$j_n(ka)=\frac{ka}{2n+1}(\frac{n}{ka}j_n(ka)+\frac{ka}{n+1}j_n(ka)),$$

this gives the condition that $$ka=n+1$$. Numerically checking this shows this is in fact wrong :(

Attempt #4

I also substituted in the expressions for $$n=1$$ into the original boundary condition to see if I noticed anything, and got this result:

$$\tan(ka)=\frac{ka}{1-(ka)^2},$$

but I think this may be a transcendental equation.

I have looked at Spherical Bessel Zeros, where they say there is no closed form for the zeros of bessel functions, however, I am wondering if there exists a closed form solution for the addition of two Bessel function. I have seen the approximation formula on Wikipedia as well as in references like Watson's Treatise on the theory of Bessel Functions, but I am trying to find a closed form solution, so I don't have to rely on numerical methods. As I recall, Bessel functions of different order do not have any overlapping zeros, as such the solution where $$j_n(ka)=0$$ and $$j_{n-1}(ka)=0$$ (the different order Bessel's share a root) does not exist.

• A similar question was asked here, but it looks like the boundary condition wasn't solved there.
Sep 22, 2020 at 4:16
• I think that I have something. Look at my answer. Sep 22, 2020 at 10:22
• I think that we could do better. By the way, Welcome to the site ! Sep 22, 2020 at 10:41
• The ratios $j_n(x)/j_{n-1}(x)$ admit an expansion in continued fractions, see <a href="dlmf.nist.gov/10.10">DLMF</a>, using $j_n=J_{n+1/2}$ because the other factors cancel. So for large $ka$ one may truncate this expansion and gets a polynomial equation for $ka$ which may be solved with polyonmial root-finding algorithms. Sep 22, 2020 at 12:07
• I have new material which will be in an edit to my answer tomorrow morning (dinner time here). Sep 22, 2020 at 15:43

I was interested by the transcendental equation. Letting $$x=ka$$, for stability reasons, I prefered to change its writng and consider that we look for the zero of function $$f(x)=\left(1-x^2\right) \sin (x)-x \cos (x)$$ Since it is odd, discarding the trivial solution $$x=0$$, focus on the positive roots which are closer and closer to multiples of $$\pi$$.

Written as series around $$x=k \pi$$, $$f(x)=-\pi k-\pi ^2 k^2 t-\frac{3\pi k}{2} t^2+\frac{\pi ^2 k^2-4}{6} t^3+\frac{7\pi k}{24} t^4+\frac{16-\pi ^2 k^2}{120} t^5-\frac{11\pi k}{720} t^6+O\left(t^{7}\right)$$ where $$t=x-k \pi$$. Now, using series reversion $$x_{(k)}=\pi k-\frac{1}{\pi k}-\frac{5}{3 \pi ^3 k^3}-\frac{73}{15 \pi ^5 k^5}-\frac{1826}{105 \pi ^7 k^7}+O\left(\frac{1}{k^9}\right)$$ which seems to work quite well.

$$\left( \begin{array}{ccc} k & \text{estimation} & \text{solution} \\ 1 & 2.747869217 & 2.743707270 \\ 2 & 6.116769339 & 6.116764264 \\ 3 & 9.316615752 & 9.316615629 \\ 4 & 12.48593738 & 12.48593737 \\ 5 & 15.64386611 & 15.64386611 \end{array} \right)$$

Edit

Considering that we look for the first non trivial zero of function $$f(x)=j_n(x)-\frac{x }{n}\,j_{n-1}(x)$$ we have $$f'(x)=\sqrt{\frac{\pi }{2}}\,\,\,\frac{x^2-n(n+1) }{n \,x^{3/2}}J_{n+\frac{1}{2}}(x)$$ $$f''(x)=\sqrt{\frac{\pi }{2}}\,\,\,\frac{(n+1) \left(x^2-n(n-1) \right)J_{n+\frac{1}{2}}(x)-x \left(x^2-n (n+1)\right)J_{n+\frac{3}{2}}(x) }{n \,x^{5/2} }$$ The first derivative cancels at $$x_*=\sqrt{n(n+1)}$$ To get an estimate, expand as series to get $$x_0=x_*+\sqrt{-2\frac {f(x_*)}{f''(x_*)}}$$ that is to say $$x_0=\sqrt{n(n+1)} +\sqrt{\sqrt{n(n+1)}\,\,\frac{ j_{n-1}\left(\sqrt{n(n+1)} \right)}{j_n\left(\sqrt{n(n+1)} \right)}-n}$$ Iterations of Newton method $$x_{k+1}=x_k+\frac{x_k \left(x_k J_{n-\frac{1}{2}}(x_k)-n J_{n+\frac{1}{2}}(x_k)\right)}{\left(x_k^2-n(n+1)\right) J_{n+\frac{1}{2}}(x_k)}$$

A few values for comparison $$\left( \begin{array}{cccc} n & x_0 & x_1 & x_2 & \text{solution} \\ 1 & 2.6691616 & 2.7445999 & 2.7437074 & 2.7437074 \\ 2 & 3.8798178 & 3.8702541 & 3.8702386 & 3.8702386 \\ 3 & 5.0331291 & 4.9740292 & 4.9734204 & 4.9734204 \\ 4 & 6.1579411 & 6.0635123 & 6.0619498 & 6.0619498 \\ 5 & 7.2649512 & 7.1428326 & 7.1402287 & 7.1402288 \\ 6 & 8.3595285 & 8.2144886 & 8.2108446 & 8.2108444 \\ 7 & 9.4448169 & 9.2801189 & 9.2754678 & 9.2754680 \\ 8 & 10.522842 & 10.340863 & 10.335248 & 10.335242 \\ 9 & 11.594999 & 11.397549 & 11.391017 & 11.391017 \\ 10 & 12.662292 & 12.450800 & 12.443395 & 12.443395 \\ 20 & 23.176232 & 22.866090 & 22.851803 & 22.851767 \\ 30 & 33.541865 & 33.167335 & 33.148172 & 33.148111 \\ 40 & 43.834852 & 43.410759 & 43.387766 & 43.387766 \\ 50 & 54.083230 & 53.618117 & 53.591940 & 53.592598 \\ 60 & 64.300908 & 63.800426 & 63.771501 & 63.771384 \\ 70 & 74.495928 & 73.964130 & 73.932777 & 73.932645 \\ 80 & 84.673409 & 84.113369 & 84.079832 & 84.079686 \\ 90 & 94.836838 & 94.250979 & 94.215451 & 94.215292 \\ 100 & 104.98871 & 104.37900 & 104.34163 & 104.34031 \\ 200 & 206.13568 & 205.34898 & 205.29815 & 205.29790 \\ 300 & 306.94240 & 306.03338 & 305.97340 & 305.97309 \\ 400 & 407.58540 & 406.57967 & 406.51253 & 406.51216 \\ 500 & 508.12878 & 507.04171 & 506.96858 & 506.96818 \\ 600 & 608.60398 & 607.44601 & 607.36768 & 607.36725 \\ 700 & 709.02905 & 707.80782 & 707.72487 & 707.72440 \\ 800 & 809.41545 & 808.13681 & 808.04967 & 807.51428 \\ 900 & 909.77091 & 908.43955 & 908.34858 & 907.81456 \\ 1000 & 1010.1010 & 1008.7207 & 1008.6262 & 1008.0934 \end{array} \right)$$

It has been numerically checked that for any $$n$$ (at least for $$n \leq 1000$$) $$f(x_0) \times f''(x_0) > 0$$ which implies, by Darboux theorem, that the solution will always be reached without any overshoot. However, since this product decreases quite fast when $$n$$ increases, we can expect more iterations for large $$n$$ (this is quite clear from the table).

I was hoping that, starting from $$x_0$$, Halley or Housholder methods would be leading to better $$x_1$$; this is not the case, the improvement being systematically quite marginal.

However, there is a way to avoid iterations (this means "almost" exact solution). Build around $$x_0$$ the $$[1,m]$$ Padé approximant which will write $$f(x) \sim \frac {f(x_0)+a^{(m)}(x-x_0) } {1+\sum_{i=1}^m b_i^{(m)} (x-x_0)^i}$$ which is at least $$O\big((x-x_0)^{m+1} \big)$$; all required coefficients are defined by the function and derivatives values at $$x=x_0$$. This gives $$x_{(m)}=x_0-\frac{f(x_0) } {a^{(m)} }$$ This has been done for this problem with $$m=1,2,3$$ and the results are better and better.

• Interesting work, I will have to compare with Python's fsolve, to see how the precision compares. I am wondering if a closed form solution even exists for variable $n$.