Finding the first N roots of transcendental equation

I need to find the first $n$ roots of the transcendental equation

\begin{equation} F(k) = J_m'(kr)Y_m'(k)-J'_m(k)Y'_m(kr) \end{equation}

for integer values of $m$ and any $r \in [0,1)$ where $J'$ and $Y'$ are derivatives of Bessel function of first and second kind. This is a (standard?) problem encountered in eigenfunction expansions for annular cylinders. Using the recurrence properties of the Bessel functions we can rewrite this equation without the derivatives, as shown in the code below.

To find the roots I'm using standard root finding algorithm fsolve in python which requires initial guesses. I have gathered that a decent good guess for the first root is going to be $m$, the order of the Bessel function we want to find roots for. From the plot I also gather that the roots are somewhat evenly spaced so that after I find the second root I can find the other guesses.

The spacing of the roots is nonlinear in $r$. Im sure there is a robust method out there that can help with the root guessing, or maybe there is a better way to go about doing this.

from scipy.special import jn, yn
from scipy.optimize import fsolve
import matplotlib.pyplot as plt
import numpy as np
from scipy.signal import argrelextrema

# Recurrence relations to relate derivative of bessels to normal bessel
Jp = lambda m,x : 0.5*(jn(m-1,x)-jn(m+1,x))
Yp = lambda m,x : 0.5*(yn(m-1,x)-yn(m+1,x))

# r is the ratio of inner/ outer radii and is [0,1)
F = lambda k,m,r : Jp(m,k*r)*Yp(m,k)-Jp(m,k)*Yp(m,k*r)

plt.figure()
k_array = np.linspace(0.1,80,5000)
m = 5

r = 0.9

F_array = F(k_array,m,r)

guesses = np.asarray([m,m/(1-r)])

roots = fsolve(F,guesses,args=(m,r,))
#iroots = argrelextrema(F_array**2,np.less)
#print iroots

plt.plot(k_array,F_array)
#plt.plot(k_array[iroots],np.zeros(len(iroots)),"ro")
plt.plot(roots,np.zeros(len(roots)),"ro")
• How big can $n$ be? Can it ever be very large, $\gtrsim 10^6$? – Kirill Apr 12 '18 at 22:35
• @Kirill No not at all, I wouldn't imagine I would need no more than 100 or so, at least of that order. – Dipole Apr 13 '18 at 0:45

I will complement @Richard Zhang 's answer (+1) with a python implementation of his suggested approach. The MATLAB package Chebfun has been partially ported in python. Actually there are two versions available: chebpy and pychebfun.

Here's an implementation of the root finding procedure with pychebfun (the approach is similar with chebpy)

from scipy.special import jn, yn
import pychebfun

# Recurrence relations to relate derivative of bessels to normal bessel
Jp = lambda m,x : 0.5*(jn(m-1,x)-jn(m+1,x))
Yp = lambda m,x : 0.5*(yn(m-1,x)-yn(m+1,x))

# r is the ratio of inner/ outer radii and is [0,1)
F = lambda k,m,r : Jp(m,k*r)*Yp(m,k)-Jp(m,k)*Yp(m,k*r)

m, r = 5, 0.9

# obtain the chebyshev approximation of the function for a specific domain
f_cheb = pychebfun.Chebfun.from_function(lambda x: F(x, m, r), domain = (0.1,80))

# obtain the roots of the function in this domain
roots = f_cheb.roots()

#plot function and roots
k_array = np.linspace(0.1,80,5000)
plt.plot(k_array, F(k_array, m, r), label = '$F$')
plt.plot(f_cheb.roots(), F(f_cheb.roots(), m, r), 'o', label = 'roots')
plt.ylim(ymin=-.1, ymax = 0.1)
plt.legend()
plt.grid()
plt.xlabel('$k$'); • Didn't realize that there's a py version of chebfun! Does it do 2D as well? – Richard Zhang Apr 12 '18 at 20:48
• @Richard Zhang Unfortunately not. – Stelios Apr 12 '18 at 20:49
• Great! (+1) I saw these two version in python and was wondering if you had an idea of which is more robust/faster? – Dipole Apr 13 '18 at 0:48
• @Jack I have only checked for "simple" tasks like this one, and they both appear the same to me. Their interface is also very similar. – Stelios Apr 13 '18 at 6:18
• @Jack, you can increase the number of polynomials used in your expansion. – nicoguaro Apr 17 '18 at 13:16

This is a root-finding problem for an analytic function over a single dimension. One standard technique is to approximate $F(k)$ as a polynomial $F(k) \approx c_0 + c_1 k + c_2 k^2 + \cdots$ using a Chebyshev approximation, and to compute the roots of the polynomial semi-analytically, e.g. by setting up the companion matrix and computing its eigenvalues. (Given that $F(k)$ is analytic, it is well-approximated by a polynomial, and the convergence of a Chebyshev polynomial approximation is exponential.) The roots found using this process can be optionally refined using Newton's method.

The above algorithm is implemented in the MATLAB package Chebfun. The reference on this algorithm is due to John P. Boyd.

• Boyd, John P. "Computing zeros on a real interval through Chebyshev expansion and polynomial rootfinding." SIAM Journal on Numerical Analysis 40.5 (2002): 1666-1682.
• Thanks! I was unaware of this method. Ill have to implement it for the present problem. – Dipole Apr 12 '18 at 20:18
• Do you need to simultaneously find the roots over $k$ and $r$? The two-dimensional root-finding problem is significantly harder. In this case, I would refer you to some of Alex Townsend's work on this problem. – Richard Zhang Apr 12 '18 at 20:37
• r is a fixed parameter in my case, so I'm guessing not - is that what you mean? – Dipole Apr 13 '18 at 0:49
• FWIW: Boyd has expanded that article to a full monograph. – J. M. Apr 13 '18 at 5:44

Complementing the answers, I would like to say that directly discretizing the differential equation might work really well. I used that approach in a previous question. I used Finite Differences and used an eigensolver to find the roots. That, besides a perturbation analysis followed by Newton's method.

The suggestion there was to use Chebyshev approximation, as well.