Finding the first N roots of transcendental equation

Question

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)[0]
#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")

Answer 1

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$');

Answer 2

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.

Answer 3

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.