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,

Jm(ωk,mr)cos(mθ)cos(cωk,mt),

where (r,θ) are polar coordinates, 0r1, t is time, m is a positive integer, c is a material parameter, and Jm is a Bessel function of the first kind. The quantity ωk,m is a resonant frequency and is a positive root of the equation

Jm(ω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
Copy to clipboard
J3(x) = besselj(3,x)
plot(J3,0,20,
    grid=:xy, legend=:none,
    xaxis=("\$x\$"), yaxis=("\$J_3(x)\$"), title="Bessel function")
Copy to clipboard

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
Copy to clipboard
pretty_table([omega J3.(omega)],["root estimate" "function value"],backend=:html)
Copy to clipboard
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")
Copy to clipboard