The Laplacian in Polar Coordinates and the Fourier-Bessel Series

We will assume a product solution of the form u_p(r,t) = W(r)G(t). Substituting this into the wave equation [@eq:wave-eq-polar-symmetric] gives $$\begin{align} \frac{\partial^2}{\partial t^2} (W G) &= c^2 \left( \frac{\partial^2}{\partial r^2} (W G) + \frac{1}{r} \frac{\partial}{\partial r} (W G) \right) \implies \\ W \frac{\partial^2 G}{\partial t^2} &= c^2 \left( G \frac{\partial^2 W}{\partial r^2} + \frac{1}{r} G \frac{\partial W}{\partial r} \right). \end{align}$$ Using Ġ = ∂G/∂t and W′ = ∂W/∂r, we rewrite this as $$ W \ddot{G} = c^2 \left( G W'' + \frac{1}{r} G W' \right), $$ which can be separated into two ordinary differential equations by dividing by c²WG: $$ \frac{\ddot{G}}{c^2 G} = \frac{W''}{W} + \frac{1}{r} \frac{W'}{W}. $$

The left-hand side depends only on t, and the right-hand side depends only on r. Therefore, both sides must be equal to a constant, which we will call  − k². This gives two ordinary differential equations, one for time and one for radius. The time equation is, letting λ ≡ ck, $$ \ddot{G} + \lambda^2 G = 0, $$ which is the familiar harmonic oscillator equation with solutions G(t) = Acos (λt) + Bsin (λt).

The radial equation is $$ W'' + \frac{1}{r} W' + k^2 W = 0, $$ which can be transformed into the Bessel equation (see example) by changing variables to s = kr. Apply the chain rule to find $$ W' = \frac{d W}{d r} = \frac{d W}{d s} \frac{d s}{d r} = k \frac{d W}{d s}, $$ and therefore W″ = k²d²W/ds². Substitute these into the radial equation to get $$ k^2 \frac{d^2 W}{d s^2} + \frac{1}{r} k \frac{d W}{d s} + k^2 W = 0 \implies \frac{d^2 W}{d s^2} + \frac{1}{s} \frac{d W}{d s} + W = 0, $$ which is the Bessel equation with ν = 0. From example, we know that the general solution to the Bessel equation is W(s) = aJ₀(s) + bY₀(s), where J₀(s) and Y₀(s) are order-0 Bessel functions of the first and second kind, respectively. The boundary condition |u| < ∞ at r = 0 requires b = 0, so the general solution for W(s) is W(s) = aJ₀(s). The boundary condition u = 0 at r = R requires J₀(kR) = 0. Zeros of the Bessel function J₀(s) are denoted by α_0, m, so the values of k are $$ k_m = \frac{\alpha_{0,m}}{R} $$ and the eigenfunctions are u_m(r,t) = (A_mcos(λ_mt)+B_msin(λ_mt))J₀(k_mr), where the eigenvalues are λ_m = ck_m.

Applying the initial condition u′(r,0) = 0 gives $$\begin{align} \left.\frac{\partial u}{\partial t}\right|_{t=0} &= -A_m \lambda_m \sin\left(\lambda_m 0\right) + B_m \lambda_m \cos\left(\lambda_m 0\right) J_0(k_m r) \implies \\ 0 &= B_m \lambda_m J_0(k_m r) \implies \\ B_m &= 0. \end{align}$$ Applying the initial condition u(r,0) = f(r) gives $$\begin{align} u(r,0) &= A_m \cos\left(\lambda_m 0\right) J_0(k_m r) \implies \\ &= A_m J_0(k_m r). \end{align}$$ We can choose A_m such that u(r,0) = f(r) in the rare cases where f(r) is a Bessel function. Instead, we will use the orthogonality of the Bessel functions to find a series representation of f(r) in terms of the eigenfunctions u_m(r,t); that is, we will find A_m such that $$ f(r) = \sum_{m=1}^\infty A_m J_0(k_m r). $$ So A_m are the Fourier-Bessel coefficients of f(r), $$ A_m = \frac{\langle f, f_m \rangle_{w=r}}{\langle f_m, f_m \rangle_{w=r}}, $$ with weight function w = r, interval [0,R], and f_m(r) = J₀(k_mr). The inner products are given from the definition of the inner product, $$\begin{align} \langle f, f_m \rangle_{w=r} &= \int_0^R f(r) f_m(r) r \, dr, \\ \langle f_m, f_m \rangle_{w=r} &= \int_0^R f_m(r)^2 r \, dr. \end{align}$$

We proceed in Python. First, load packages.

import numpy as np
import sympy as sp
import scipy
from scipy.special import jn_zeros
import matplotlib.pyplot as plt

The inner products are computed as follows:

r = sp.symbols("r", real=True, nonnegative=True)
k_m, R = sp.symbols("k_m, R", real=True, positive=True)
f = 1 - (r/R)**3  # Initial condition
f_m = sp.besselj(0, k_m * r)  # Eigenfunction
w = r  # Weight function of the inner product
inner_f_fm = sp.integrate(f * f_m * w, (r, 0, R))  # $\langle f, f_m \rangle$
inner_fm_fm = sp.integrate(f_m**2 * w, (r, 0, R))  # $\langle f_m, f_m \rangle$
A_m = (inner_f_fm / inner_fm_fm).simplify()

The values of k_m, eigenvalues λ_m = ck_m, and the series components A_m can be computed from the zeros of the Bessel function zeros α_0, m, where m = 1, 2, 3, … as follows:

c = sp.symbols("c", real=True, positive=True)
n_zeros = 5  # Number of Bessel function zeros
params = {R: 1, c: 2*np.pi}

The partial sum of the series can be defined as the sum of the first n terms of the series:

def f_partial_sum_t(r, t, n, A_m, params):
    alphas = jn_zeros(0, n)  # Bessel function zeros
    ks = alphas / params[R]
    lambdas = params[c] * ks
    A_m_ = np.zeros(n)  # Numerical vals of the series comp.
    f_partial = np.zeros((len(r), len(t)))
    for i in range(n):
        A_m_[i] = A_m.subs(k_m, ks[i]).subs(params).evalf()
        f_m = scipy.special.jn(0, ks[i] * r[:, None])
        f_partial += A_m_[i] * f_m * np.cos(lambdas[i] * t)
    return f_partial

Plot the partial sum of the series and compare it to the initial condition.

r_plt = np.linspace(-params[R], params[R], 101)
t_plt = np.array([0])
f_num = 1 - (np.abs(r_plt)/params[R])**3
f_partial = f_partial_sum_t(r_plt, t_plt, n_zeros, A_m, params)
fig, ax = plt.subplots()
ax.plot(r_plt, f_num, label="$f$")
ax.plot(r_plt, f_partial, label="$f_{\mathrm{partial}}$")
ax.set_xlabel('$r$')
ax.legend()
plt.draw()

The result is quite good, even with only five terms in the series.

We can now plot the partial sum over one period.

cmap = plt.get_cmap('viridis')
n_snapshots = 101
t_plt = np.linspace(0, 1, n_snapshots)
r_plt = np.linspace(-params[R], params[R], 201)
f_partial = f_partial_sum_t(r_plt, t_plt, n_zeros, A_m, params)
fig, ax = plt.subplots()
for i in range(n_snapshots):
    # Label 5 snapshots
    label = None
    if i % 20 == 0:
        label = f"$t = {t_plt[i]:.2f}$"
    ax.plot(r_plt, f_partial[:, i], label=label,
            color=cmap(i / n_snapshots), alpha=1)
ax.set_xlabel('$r$')
ax.set_ylabel('Displacement u(r, t)')
ax.legend(loc='upper right')
plt.show()