Roots of Bessel functions

In the theory of vibrations of a circular drum, the displacement of the drumhead can be expressed in terms of pure harmonic modes,

\[J_m(\omega_{k,m} r) \cos(m\theta) \cos(c \omega_{k,m} t),\]

where \((r,\theta)\) are polar coordinates, \(0\le r\le 1\), \(t\) is time, \(m\) is a positive integer, \(c\) is a material parameter, and \(J_m\) is a Bessel function of the first kind. The quantity \(\omega_{k,m}\) is a resonant frequency and is a positive root of the equation

\[J_m(\omega_{k,m}) = 0,\]

which states that the drumhead is clamped around the rim. Tabulating approximations to the zeros of Bessel functions has occupied countless mathematician-hours throughout the centuries.

using FundamentalsNumericalComputation,SpecialFunctions

J3(x) = besselj(3,x)
plot(J3,0,20,
    grid=:xy, legend=:none,
    xaxis=("\$x\$"), yaxis=("\$J_3(x)\$"), title="Bessel function")

From the graph we see roots near 6, 10, 13, 16, and 19. We use nlsolve from the NLsolve package to find these roots accurately. (It uses vector variables, so we have to adapt it for use with scalars.)

omega = []
for guess = [6.,10.,13.,16.,19.]
    s = nlsolve(x->besselj(3,x[1]),[guess])
    omega = [omega;s.zero]
end

pretty_table([omega J3.(omega)],["root estimate" "function value"],backend=:html)
root estimate function value
6.380161895923975 2.4702462297909733e-15
9.761023129981334 -8.375244942016025e-14
13.015200721696122 5.048461648726743e-13
16.223466160318768 6.938893903907228e-18
19.40941522643161 6.126904539272005e-13 
scatter!(omega,J3.(omega),title="Bessel function with roots")
Chapter 4 (Rootfinding) Conditioning of roots