Boundary Value Problems

This is a sturm-liouville problem, so we know the eigenvalues are real. The well-known general solution to the ODE is $$X(x) = \begin{cases} k_1 + k_2 x & \lambda = 0 \\ k_1 e^{j\sqrt{\lambda} x} + k_2 e^{-j\sqrt{\lambda} x} & \text{otherwise} \end{cases}$$ with real constants k₁, k₂. The solution must also satisfy the boundary conditions. Let’s apply them to the case of λ = 0 first: $$\begin{aligned} X(0) &= 0 \Rightarrow k_1 + k_2 (0) = 0 \Rightarrow k_1 = 0 \\ X(L) &= 0 \Rightarrow k_1 + k_2 (L) = 0 \Rightarrow k_2 = -k_1/L. \end{aligned}$$ Together, these imply k₁ = k₂ = 0, which gives the trivial solution X(x) = 0, in which we aren’t interested. We say, then, that for nontrivial solutions, λ ≠ 0. Now let’s check λ < 0. The solution becomes $$\begin{aligned} X(x) &= k_1 e^{-\sqrt{\abs{\lambda}} x} + k_2 e^{\sqrt{\abs{\lambda}} x} \\ &= k_3 \cosh(\sqrt{\abs{\lambda}}x) + k_4 \sinh(\sqrt{\abs{\lambda}}x) \end{aligned}$$ where k₃ and k₄ are real constants. Again applying the boundary conditions: $$\begin{aligned} X(0) &= 0 \Rightarrow k_3 \cosh(0) + k_4 \sinh(0) = 0 \Rightarrow k_3 + 0 = 0 \Rightarrow k_3 = 0 \\ X(L) &= 0 \Rightarrow 0 \cosh(\sqrt{\abs{\lambda}}L) + k_4 \sinh(\sqrt{\abs{\lambda}}L) = 0 \Rightarrow k_4 \sinh(\sqrt{\abs{\lambda}}L) = 0. \end{aligned}$$ However, $\sinh(\sqrt{\abs{\lambda}}L) \ne 0$ for L > 0, so k₄ = k₃ = 0—again, the trivial solution. Now let’s try λ > 0. The solution can be written $$\begin{aligned} X(x) = k_5 \cos(\sqrt{\lambda}x) + k_6 \sin(\sqrt{\lambda}x). \end{aligned}$$ Applying the boundary conditions for this case: $$\begin{aligned} X(0) &= 0 \Rightarrow k_5 \cos(0) + k_6 \sin(0) = 0 \Rightarrow k_5 + 0 = 0 \Rightarrow k_5 = 0 \\ X(L) &= 0 \Rightarrow 0 \cos(\sqrt{\lambda}L) + k_6 \sin(\sqrt{\lambda}L) = 0 \Rightarrow k_6 \sin(\sqrt{\lambda}L) = 0. \end{aligned}$$ Now, $\sin(\sqrt{\lambda}L) = 0$ for $$\begin{align} \sqrt{\lambda}L &= n \pi \Rightarrow \nonumber \\ \lambda &= \left(\frac{n \pi} {L}\right)^2. \tag{$n \in \mathbb{Z_+}$} \end{align}$$ Therefore, the only nontrivial solutions that satisfy both the ODE and the boundary conditions are the eigenfunctions $$\begin{align} X_n(x) &= \sin\left(\sqrt{\lambda_n} x\right) \\ &= \sin\left(\frac{n \pi} {L} x\right) \end{align}$$ with corresponding eigenvalues $$\begin{aligned} \lambda_n &= \left(\frac{n \pi} {L}\right)^2. \end{aligned}$$ Note that because λ > 0, λ₁ is the lowest eigenvalue.

Plotting the eigenfunctions

import numpy as np
import matplotlib.pyplot as plt

Set L = 1 and compute values for the first four eigenvalues lambda_n and eigenfunctions X_n.

L = 1
x = np.linspace(0, L, 100)
n = np.linspace(1, 4, 4, dtype=int)
lambda_n = (n*np.pi/L)**2
X_n = np.zeros([len(n), len(x)])
for i,n_i in enumerate(n):
  X_n[i, :] = np.sin(np.sqrt(lambda_n[i])*x)

Plot the eigenfunctions.

fig, ax = plt.subplots()
for i, n_i in enumerate(n):
  ax.plot(x, X_n[i,:], linewidth=2,label='$n = '+str(n_i)+'$')
plt.legend()
plt.show()

We see that the fourth of the S-L theorems appears true: n − 1 zeros of X_n exist on the open interval (0,1).

We proceed to show that the Bessel equation with the given boundary conditions is an irregular Sturm-Liouville problem by first showing that the Bessel equation can be written in the form of a Sturm-Liouville problem ODE, then showing that the boundary conditions are not regular. Dividing the Bessel equation by s gives $$\begin{align} s y'' + y' + \frac{s^2 - \nu^2}{s} y &= 0 \implies \tag{z \ne 0} \\ \frac{d}{d s} \left( s y' \right) + \left(s - \nu^2/s\right) y &= 0. \end{align}$$ So we see that this is equivalent to the Sturm-Liouville problem’s ODE, $$\begin{align} \frac{d} {d x} \left(p X'\right) + \left(\lambda \sigma + q\right) X = 0, \end{align}$$ with x = s, X = y, p(s) = s, q(s) =  − ν²/s, σ(s) = s, and λ = 1.

Now let us consider the boundary conditions. The second boundary condition cannot be written as a linear combination of y and y′ at s = 0, therefore this is an irregular Sturm-Liouville problem.

Solutions for Bessel’s equation are Bessel functions of the first kind, J_ν(s) (section 6.4), and Bessel functions of the second kind, Y_ν(s) [@kreyszig2011,§ 5.5]. For ν = 0, the Bessel equation simplifies to s²y″ + sy′ + s²y = 0 and the solutions are zeroth-order Bessel functions of the first kind, J₀(s), and of the second kind, Y₀(s). That is, the general solution is y(s) = aJ₀(s) + bY₀(s). However, Y₀(s) is singular at s = 0, so the boundary condition |y(0)| < ∞ requires b = 0. Therefore, the eigenfunctions are y_n(s) = J₀(s) = J₀(k_nr) for eigenvalues λ_n = k_n² and n ∈ ℤ₊. On the boundary, y_n(kR) = 0 implies that k_nR are the zeros of J₀(s), what we called α_0, n in example. Therefore, the eigenvalues are $$ \lambda_n = \left(\frac{\alpha_{0,n}}{R}\right)^2. $$

Plotting the eigenfunctions

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 eigenvalues can be computed from the zeros of the Bessel function zeros α_0, n, where n = 1, 2, 3, … as follows:

n = sp.symbols('n', integer=True, positive=True)
k, R = sp.symbols('k, R', positive=True, real=True)
r = sp.symbols('r', nonnegative=True)
J_0 = sp.besselj(0, k * r)
N_zeros = 5
alpha_0_n = jn_zeros(0, N_zeros)
params = {R: 1}  # Set R = 1
lambda_n_ = (alpha_0_n / params[R])**2
print(f"The first {N_zeros} eigenvalues are:")
print(lambda_n_)

The first 5 eigenvalues are:
[  5.78318596  30.47126234  74.88700679 139.04028443 222.93230362]

The eigenfunctions are given by the Bessel functions of the first kind J₀(k_nr).

def k_n(lambda_n): return np.sqrt(lambda_n)
k_n_ = k_n(lambda_n_)

Plot the eigenfunctions.

r_plt = np.linspace(0, 1, 101)
print(J_0)
y_n_fun = sp.lambdify((r, k), J_0, modules=['numpy', 'scipy'])
fig, ax = plt.subplots()
for i in range(5):
    k_i = k_n_[i]
    ax.plot(r_plt, y_n_fun(r_plt, k_i), label=f'$y_{i+1}(k_{i+1} r)$')
ax.spines['bottom'].set_position('zero')
ax.spines['left'].set_position('zero')
ax.set_xlabel('$r$')
ax.set_ylabel('Eigenfunctions $y_n(k_n r)$')
ax.legend()
plt.show()

The plot shows the eigenfunctions y_n(k_nr) for the first few eigenvalues. Note that the boundary conditions are satisfied and that the eigenfunctions appear to be radial modes of vibration for a circular membrane.