Software 11.6.1 (The drum’s frequencies). The zeros \(j_{mn}\) of the Bessel functions for \(m \leq 5\), \(n \leq 3\), sorted, and the frequencies of the first eight modes relative to the fundamental:

>>> import numpy as np
>>> from scipy.special import jn_zeros, jv
>>> modes = [(m, n, jn_zeros(m, n)[-1])
...          for m in range(6) for n in range(1, 4)]
>>> modes.sort(key=lambda t: t[2])
>>> base = modes[0][2]                     # j_01 = 2.4048: the fundamental
>>> for m, n, j in modes[:8]:  # frequency ratios of a circular drum
...     print((m, n), round(j, 4), round(j / base, 3))
...
(0, 1) 2.4048 1.0
(1, 1) 3.8317 1.593
(2, 1) 5.1356 2.136
(0, 2) 5.5201 2.295
(3, 1) 6.3802 2.653
(1, 2) 7.0156 2.917
(4, 1) 7.5883 3.155
(2, 2) 8.4172 3.5

Reading the output. The overtones of a drum are at \(1.593\), \(2.136\), \(2.295\), \(2.653\), \(2.917\), \(3.155\) and \(3.5\) times the fundamental, irregularly spaced. A string’s are at \(2, 3, 4, \ldots \).

MATLAB/Octave: fzero(@(x) besselj(m,x), guess); there is no built-in for the zeros.

Software 11.6.2 (Cooling a cylinder). The Bessel series for a cylinder of unit radius with \(\kappa = 1\) and \(u_0 = 1\), summed to \(60\) terms: first a check of orthogonality, then the centre temperature at three times against the first term alone:

>>> from scipy.integrate import quad
>>> j0 = jn_zeros(0, 60)        # cooling cylinder: a = kappa = 1, u0 = 1
>>> c = 2 / (j0 * jv(1, j0))               # coefficients of 1 in J_0(j_n r)
>>> quad(lambda r: r * jv(0, j0[0]*r) * jv(0, j0[1]*r), 0, 1)[0]  # orthogonal
5.204170427930421e-18
>>> def centre(t, N=60):
...     return np.sum(c[:N] * np.exp(-j0[:N]**2 * t))
...
>>> for t in (0.01, 0.1, 0.5):
...     print(t, round(centre(t), 8), round(c[0] * np.exp(-j0[0]**2 * t), 8))
...
0.01 1.0 1.51195754
0.1 0.84835511 0.89845241
0.5 0.08888972 0.08888997

Reading the output. The integral of \(J_0(j_{01}r)J_0(j_{02}r)r\) is zero to rounding. At \(t = 0.01\) the centre is still at \(1\), while the one-term value is \(1.51\); at \(t = 0.1\) the two differ by \(6\%\); by \(t = 0.5\) they agree to seven digits, the temperature having fallen to \(0.089\).

MATLAB/Octave: sum(c.*exp(-j0.^2*t)) with c = 2./(j0.*besselj(1,j0)).

Software 11.6.3 (Legendre series and Kelvin’s estimate). The coefficients \(A_l\) of Example 11.4.1 for \(l = 1, 3, \ldots , 11\), the six-term series on the axis at \(z = \frac 12\) against the closed form, and Kelvin’s age of the Earth:

>>> import sympy as sp
>>> mu = sp.symbols('mu')
>>> A = [(2*l + 1) * sp.integrate(sp.legendre(l, mu), (mu, 0, 1))
...      for l in range(1, 12, 2)]  # hemispheres at +1 and -1, odd l only
>>> A
[3/2, -7/8, 11/16, -75/128, 133/256, -483/1024]
>>> z = 0.5  # on the axis, from the Poisson kernel
>>> series = sum(float(a) * z**l for a, l in zip(A, range(1, 12, 2)))
>>> exact = (1 - (1 - z**2) / np.sqrt(1 + z**2)) / z
>>> round(series, 6), round(exact, 6)
(0.658316, 0.658359)
>>> T0, G, kappa = 3900.0, 0.0365, 1.2e-6  # Kelvin: K, K/m, m^2/s
>>> t = T0**2 / (np.pi * kappa * G**2)  # gradient T0 / sqrt(pi kappa t) = G
>>> round(t / 3.156e7 / 1e6, 1)            # in millions of years
96.0

Reading the output. The coefficients \(\frac 32, -\frac 78, \frac {11}{16}, -\frac {75}{128}, \ldots \) decay slowly, like \(l^{-1/2}\), because the boundary data jump at the equator; but at \(r = \frac 12\) the factor \(r^l\) makes six terms accurate to \(4 \times 10^{-5}\). Kelvin’s numbers give \(96.0\) million years.

MATLAB/Octave: legendre(l, mu) returns the associated functions; the first row is \(P_l\).