Ordinary Differential Equations (5): Series Solutions and Special Functions

In the world of differential equations, some equations cannot be expressed in terms of elementary functions, but this doesn't mean we're helpless. Power series solutions open a new door — by expanding solutions as infinite series, we can not only solve a broader class of equations but also give birth to special functions like Bessel functions and Legendre polynomials that appear everywhere in physics. This chapter will guide you through this exquisite mathematical landscape.

Starting with a Simple Problem

Imagine you're studying heat conduction and need to find the temperature distribution inside a cylinder. After separation of variables, you encounter this equation:

This is the famous Bessel equation. Trying all the methods we've learned before — separation of variables, integrating factors, variation of parameters — none work. The coefficients of this equation are functions of, not constants, so traditional methods fail.

Fortunately, power series solutions can save us.

The Core Idea of Power Series Solutions

Basic Principle

Core assumption: Assume the solution can be expressed as a power series

Solution steps: 1. Substitute the power series forinto the differential equation 2. Write all terms in the same power form 3. Set the coefficient of each power to zero 4. Obtain a recurrence relation for coefficients 5. Use initial conditions to determine, etc.

Review: A Familiar Example

Let's use the power series method to re-solve a simple equation:Assume, then.

Substituting into the equation:Adjusting the summation index, leton the left:Comparing coefficients of:Fromwe know, and by recurrence:Therefore:The power series method has "rediscovered" the exponential function!

Ordinary Points and Singular Points

Definitions

Consider a second-order linear equation:Converting to standard form:where,.

Definitions: - Ifandare analytic (can be expanded as convergent power series) at, thenis called an ordinary point - Otherwise,is called a singular point

Power Series Solutions at Ordinary Points

Theorem: Ifis an ordinary point, then the equation has two linearly independent power series solutions near:with radius of convergence at least equal to the distance fromto the nearest singular point.

Example: The Airy Equation

Airy equation:This is a core equation in WKB approximation in quantum mechanics and rainbow theory in optics.

Analysis:,, both are entire functions, sois an ordinary point.

Solution: Let, then.

Substituting into the equation:Adjusting indices to match powers of. For the first sum, let; for the second, let:Separating theterm:Setting each coefficient to zero: -: -:This recurrence relation means the coefficients divide into three groups: - Starting with: - Starting with: - Starting with:Computing a few terms:General solution:whereandare Airy functions— this is the birth of two special functions!

import numpy as np
import matplotlib.pyplot as plt
from scipy.special import airy

# Airy function visualization
x = np.linspace(-15, 5, 1000)
Ai, Aip, Bi, Bip = airy(x)

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 5))

# Airy functions
ax1.plot(x, Ai, 'b-', linewidth=2, label='Ai(x)')
ax1.plot(x, Bi, 'r-', linewidth=2, label='Bi(x)')
ax1.axhline(y=0, color='k', linestyle='-', linewidth=0.5)
ax1.axvline(x=0, color='k', linestyle='-', linewidth=0.5)
ax1.set_xlabel('x', fontsize=12)
ax1.set_ylabel('y', fontsize=12)
ax1.set_title('Airy Functions: Solutions of y\'\' - xy = 0', fontsize=14, fontweight='bold')
ax1.legend(fontsize=11)
ax1.grid(True, alpha=0.3)
ax1.set_ylim(-0.6, 1.2)

# Airy functions in quantum tunneling
# Linear potential V(x) = V0 * x
x_qm = np.linspace(-5, 3, 500)
V = x_qm  # Linear potential (normalized)
psi, _, _, _ = airy(x_qm)  # Wavefunction proportional to Ai(x)
psi_prob = np.abs(psi)**2  # Probability density

ax2.fill_between(x_qm, 0, V, alpha=0.3, color='gray', label='Potential V(x)')
ax2.plot(x_qm, psi_prob * 5, 'b-', linewidth=2, label='|ψ(x)|² (scaled)')
ax2.axvline(x=0, color='r', linestyle='--', linewidth=1.5, label='Classical turning point')
ax2.set_xlabel('Position x', fontsize=12)
ax2.set_ylabel('Energy / Probability', fontsize=12)
ax2.set_title('Quantum Tunneling Near Linear Potential Barrier', fontsize=14, fontweight='bold')
ax2.legend(fontsize=10)
ax2.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('airy_functions.png', dpi=150, bbox_inches='tight')
plt.show()

print("Ai(0) =", airy(0)[0])
print("Bi(0) =", airy(0)[2])

Regular Singular Points and the Frobenius Method

What is a Regular Singular Point?

Not all singular points are equally "bad." Some can be handled by an improved power series method.

Definition: Ifis a singular point, butare analytic at, thenis called a regular singular point.

Otherwise, it's called an irregular singular point.

The Frobenius Method

Idea: At a regular singular point, try a generalized power series solution:wherecan be any real number (or even complex).

Steps: 1. Substitute the generalized power series into the equation 2. Set the coefficient of the lowest power to zero to get the indicial equation 3. Solve the indicial equation for possible values of$rr$

The Indicial Equation

For an equation in standard form:Let,.

Indicial equation:Let the two roots be(with convention).

Three Cases

Depending on the nature of, there are three cases:

Case 1:is not an integer Two linearly independent solutions:

Case 2:

Case 3:is a positive integerwheremay be zero.

Bessel Equation and Bessel Functions

The Bessel EquationThis is one of the most important equations in physics, appearing in:

• Laplace equation in cylindrical coordinates
• Vibration of circular drumheads
• Electromagnetic wave propagation in cylindrical waveguides
• Heat conduction problems

Analysisis a regular singular point. Converting to standard form:So,.

At:,.

Indicial equation:Solving gives.

Bessel Functions of the First KindTaking(), let.

Substituting into the equation, after tedious but straightforward calculations, we get the recurrence relation:Usually we take,.

Result:

Bessel Functions of the Second Kind(Neumann Functions)

Whenis not an integer,is also a solution and is linearly independent from.

Whenis an integer,, so we need to construct a second solution.

Neumann function definition:Whenis an integer, it's defined by a limit.

Properties of Bessel Functions

Recurrence relations:$

$

Derivative relations:$

Asymptotic behavior: As:As:

import numpy as np
import matplotlib.pyplot as plt
from scipy.special import jv, yv

# Bessel function visualization
x = np.linspace(0.01, 20, 1000)

fig, axes = plt.subplots(2, 2, figsize=(14, 10))

# Bessel functions of the first kind
ax1 = axes[0, 0]
for n in range(5):
    ax1.plot(x, jv(n, x), linewidth=2, label=f'$J_{n}(x)$')
ax1.axhline(y=0, color='k', linestyle='-', linewidth=0.5)
ax1.set_xlabel('x', fontsize=12)
ax1.set_ylabel('$J_n(x)$', fontsize=12)
ax1.set_title('Bessel Functions of the First Kind', fontsize=14, fontweight='bold')
ax1.legend(fontsize=10)
ax1.grid(True, alpha=0.3)
ax1.set_ylim(-0.5, 1.1)

# Bessel functions of the second kind
ax2 = axes[0, 1]
for n in range(4):
    ax2.plot(x, yv(n, x), linewidth=2, label=f'$Y_{n}(x)$')
ax2.axhline(y=0, color='k', linestyle='-', linewidth=0.5)
ax2.set_xlabel('x', fontsize=12)
ax2.set_ylabel('$Y_n(x)$', fontsize=12)
ax2.set_title('Bessel Functions of the Second Kind (Neumann)', fontsize=14, fontweight='bold')
ax2.legend(fontsize=10)
ax2.grid(True, alpha=0.3)
ax2.set_ylim(-1.5, 0.6)

# Zeros of Bessel functions (key for drum vibration modes)
ax3 = axes[1, 0]
from scipy.special import jn_zeros
for n in range(4):
    zeros = jn_zeros(n, 5)
    ax3.scatter(zeros, [n]*5, s=100, label=f'$J_{n}$zeros')
    for z in zeros:
        ax3.axvline(x=z, color=f'C{n}', alpha=0.3, linestyle='--')
ax3.set_xlabel('x', fontsize=12)
ax3.set_ylabel('n (order)', fontsize=12)
ax3.set_title('Zeros of Bessel Functions (Drum Vibration Modes)', fontsize=14, fontweight='bold')
ax3.legend(fontsize=10)
ax3.grid(True, alpha=0.3)

# Circular drum vibration mode
ax4 = axes[1, 1]
theta = np.linspace(0, 2*np.pi, 100)
r = np.linspace(0, 1, 50)
R, Theta = np.meshgrid(r, theta)
X = R * np.cos(Theta)
Y = R * np.sin(Theta)

# A vibration mode (n=1, m=1)
n, m = 1, 1
zero_nm = jn_zeros(n, m)[-1]
Z = jv(n, zero_nm * R) * np.cos(n * Theta)

contour = ax4.contourf(X, Y, Z, levels=20, cmap='RdBu')
ax4.set_xlabel('x', fontsize=12)
ax4.set_ylabel('y', fontsize=12)
ax4.set_title(f'Circular Drum Vibration Mode (n={n}, m={m})', fontsize=14, fontweight='bold')
ax4.set_aspect('equal')
plt.colorbar(contour, ax=ax4)

plt.tight_layout()
plt.savefig('bessel_functions.png', dpi=150, bbox_inches='tight')
plt.show()

Legendre Equation and Legendre Polynomials

The Legendre EquationOr in standard form:This appears everywhere in spherical coordinates:

• Hydrogen atom wavefunctions
• Gravitational field expansion of the Earth
• Antenna radiation patterns

Analysisis an ordinary point,are regular singular points.

Expanding near: Let, substituting into the equation gives the recurrence relation:

Legendre Polynomials

Whenis a non-negative integer, the series truncates to a polynomial!

Because when,, and all subsequent even (or odd) coefficients are zero.

Legendre polynomials are these truncated solutions, normalized so that.

First few Legendre polynomials:$

$

$

$

Rodrigues' FormulaThis elegant formula directly gives Legendre polynomials of any order.

Orthogonality

Legendre polynomials are orthogonal onwith weight function:

Physical significance: Any "reasonable" function oncan be expanded as a Legendre seriesThis is crucial in multipole expansions!

import numpy as np
import matplotlib.pyplot as plt
from scipy.special import legendre

# Legendre polynomial visualization
x = np.linspace(-1, 1, 500)

fig, axes = plt.subplots(2, 2, figsize=(14, 10))

# First few Legendre polynomials
ax1 = axes[0, 0]
for n in range(6):
    Pn = legendre(n)
    ax1.plot(x, Pn(x), linewidth=2, label=f'$P_{n}(x)$')
ax1.axhline(y=0, color='k', linestyle='-', linewidth=0.5)
ax1.axvline(x=0, color='k', linestyle='-', linewidth=0.5)
ax1.set_xlabel('x', fontsize=12)
ax1.set_ylabel('$P_n(x)$', fontsize=12)
ax1.set_title('Legendre Polynomials', fontsize=14, fontweight='bold')
ax1.legend(fontsize=10, loc='lower right')
ax1.grid(True, alpha=0.3)

# Orthogonality verification
ax2 = axes[0, 1]
n_max = 6
ortho_matrix = np.zeros((n_max, n_max))
for i in range(n_max):
    for j in range(n_max):
        Pi = legendre(i)
        Pj = legendre(j)
        integral, _ = np.polynomial.legendre.leggauss(50)
        # Numerical integration
        from scipy.integrate import quad
        result, _ = quad(lambda t: Pi(t) * Pj(t), -1, 1)
        ortho_matrix[i, j] = result

im = ax2.imshow(ortho_matrix, cmap='coolwarm', vmin=-0.5, vmax=2.5)
ax2.set_xlabel('n', fontsize=12)
ax2.set_ylabel('m', fontsize=12)
ax2.set_title('Orthogonality Matrix$\\int_{-1}^{1} P_m P_n dx$', fontsize=14, fontweight='bold')
plt.colorbar(im, ax=ax2)
ax2.set_xticks(range(n_max))
ax2.set_yticks(range(n_max))

# Function expansion example: approximating step function with Legendre polynomials
ax3 = axes[1, 0]
def step_function(x):
    return np.where(x >= 0, 1, -1)

# Compute expansion coefficients
def legendre_coefficients(f, n_terms):
    coeffs = []
    for n in range(n_terms):
        Pn = legendre(n)
        integral, _ = quad(lambda t: f(t) * Pn(t), -1, 1)
        c_n = (2*n + 1) / 2 * integral
        coeffs.append(c_n)
    return coeffs

# Approximation
n_terms_list = [3, 7, 15, 31]
ax3.plot(x, step_function(x), 'k-', linewidth=3, label='Step function')

for n_terms in n_terms_list:
    coeffs = legendre_coefficients(step_function, n_terms)
    y_approx = np.zeros_like(x)
    for n, c in enumerate(coeffs):
        Pn = legendre(n)
        y_approx += c * Pn(x)
    ax3.plot(x, y_approx, linewidth=1.5, label=f'{n_terms} terms', alpha=0.7)

ax3.set_xlabel('x', fontsize=12)
ax3.set_ylabel('y', fontsize=12)
ax3.set_title('Approximating Step Function with Legendre Series', fontsize=14, fontweight='bold')
ax3.legend(fontsize=10)
ax3.grid(True, alpha=0.3)
ax3.set_ylim(-1.5, 1.5)

# Spherical harmonics visualization (application of Legendre polynomials)
ax4 = axes[1, 1]
theta = np.linspace(0, np.pi, 100)
phi = np.linspace(0, 2*np.pi, 200)
Theta, Phi = np.meshgrid(theta, phi)

# Spherical harmonic Y_2^0(θ, φ) ∝ P_2(cos θ)
from scipy.special import sph_harm
l, m = 2, 0
Y = sph_harm(m, l, Phi, Theta)
Y_real = np.real(Y)

# Convert to Cartesian coordinates
r = np.abs(Y_real)
X = r * np.sin(Theta) * np.cos(Phi)
Y_coord = r * np.sin(Theta) * np.sin(Phi)
Z = r * np.cos(Theta)

ax4.contourf(Phi[:, 0] * 180/np.pi, Theta[0, :] * 180/np.pi, Y_real.T, 
             levels=20, cmap='RdBu')
ax4.set_xlabel('φ (degrees)', fontsize=12)
ax4.set_ylabel('θ (degrees)', fontsize=12)
ax4.set_title(f'Spherical Harmonic$Y_{l}^{m}$', fontsize=14, fontweight='bold')

plt.tight_layout()
plt.savefig('legendre_polynomials.png', dpi=150, bbox_inches='tight')
plt.show()

Hermite Equation and Hermite Polynomials

The Quantum Harmonic Oscillator

In quantum mechanics, the one-dimensional harmonic oscillator's Schr ö dinger equation simplifies to the Hermite equation:Whenis a non-negative integer, there are polynomial solutions — the Hermite polynomials.

Definition of Hermite Polynomials

Rodrigues' formula:First few Hermite polynomials:$

$

$

$

Orthogonality

Hermite polynomials are orthogonal with weight functionon:

Quantum Harmonic Oscillator Wavefunctions

Energy eigenstates of the harmonic oscillator:Energy levels:This is the famous energy quantization!

import numpy as np
import matplotlib.pyplot as plt
from scipy.special import hermite
from math import factorial

# Hermite polynomials and quantum harmonic oscillator
fig, axes = plt.subplots(2, 2, figsize=(14, 10))

x = np.linspace(-4, 4, 500)

# Hermite polynomials
ax1 = axes[0, 0]
for n in range(5):
    Hn = hermite(n)
    ax1.plot(x, Hn(x), linewidth=2, label=f'$H_{n}(x)$')
ax1.axhline(y=0, color='k', linestyle='-', linewidth=0.5)
ax1.axvline(x=0, color='k', linestyle='-', linewidth=0.5)
ax1.set_xlabel('x', fontsize=12)
ax1.set_ylabel('$H_n(x)$', fontsize=12)
ax1.set_title('Hermite Polynomials', fontsize=14, fontweight='bold')
ax1.legend(fontsize=10)
ax1.grid(True, alpha=0.3)
ax1.set_ylim(-50, 50)

# Quantum harmonic oscillator wavefunctions
ax2 = axes[0, 1]
# Setting ℏ = m = ω = 1 (natural units)
def psi_n(n, x):
    Hn = hermite(n)
    norm = 1 / np.sqrt(2**n * factorial(n)) * (1/np.pi)**0.25
    return norm * Hn(x) * np.exp(-x**2 / 2)

for n in range(5):
    psi = psi_n(n, x)
    # Offset along y-axis for visualization
    ax2.plot(x, psi + n, linewidth=2, label=f'n={n},$E_n$={n+0.5}')
    ax2.fill_between(x, n, psi + n, alpha=0.3)
    ax2.axhline(y=n, color='gray', linestyle='--', linewidth=0.5)

ax2.set_xlabel('Position x', fontsize=12)
ax2.set_ylabel('$\\psi_n(x)$+ Energy offset', fontsize=12)
ax2.set_title('Quantum Harmonic Oscillator Wavefunctions', fontsize=14, fontweight='bold')
ax2.grid(True, alpha=0.3)

# Probability densities
ax3 = axes[1, 0]
for n in range(4):
    psi = psi_n(n, x)
    prob = np.abs(psi)**2
    ax3.plot(x, prob, linewidth=2, label=f'n={n}')

# Add classical probability distribution comparison (approximation for n=10)
n_classical = 10
A = np.sqrt(2*n_classical + 1)  # Classical amplitude
x_classical = np.linspace(-A, A, 500)
prob_classical = 1 / (np.pi * np.sqrt(A**2 - x_classical**2 + 0.01))
ax3.plot(x_classical, prob_classical / 10, 'k--', linewidth=2, 
         label='Classical (scaled)', alpha=0.7)

ax3.set_xlabel('Position x', fontsize=12)
ax3.set_ylabel('$|\\psi_n(x)|^2$', fontsize=12)
ax3.set_title('Probability Densities', fontsize=14, fontweight='bold')
ax3.legend(fontsize=10)
ax3.grid(True, alpha=0.3)

# Energy level diagram
ax4 = axes[1, 1]
n_levels = 8
for n in range(n_levels):
    E = n + 0.5
    ax4.hlines(E, 0.2, 0.8, linewidth=3, color=f'C{n % 10}')
    ax4.text(0.85, E, f'n={n},$E_{n}$={E}', fontsize=10, va='center')

# Draw potential curve
x_pot = np.linspace(-3, 3, 100)
V = 0.5 * x_pot**2
ax4.plot(x_pot + 0.5, V, 'k-', linewidth=2, label='$V(x) = \\frac{1}{2}x^2$')

ax4.set_xlim(-0.5, 1.5)
ax4.set_ylim(-0.5, 8)
ax4.set_xlabel('Position / Level', fontsize=12)
ax4.set_ylabel('Energy', fontsize=12)
ax4.set_title('Energy Levels of Quantum Harmonic Oscillator', fontsize=14, fontweight='bold')
ax4.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('hermite_quantum.png', dpi=150, bbox_inches='tight')
plt.show()

Practical Application: Heat Conduction and Fourier-Bessel Series

Problem Setup

Consider a cylindrical metal rod of radius, with initial temperature distribution, and surface maintained at zero temperature. Find how the temperature changes with time.

Mathematical Model

Heat conduction equation (cylindrical coordinates, assuming temperature depends only onand):Boundary condition:Initial condition:

Separation of Variables

Let, substituting gives:Time part:Spatial part:Let, this is exactly the zeroth-order Bessel equation!

Solution

General solution:Sincediverges at, take.

Boundary conditionrequires.

Let the zeros ofbe, then.

Final solution:whereare determined by initial conditions and orthogonality of Bessel functions.

import numpy as np
import matplotlib.pyplot as plt
from scipy.special import jv, jn_zeros
from scipy.integrate import quad

# Cylindrical heat conduction problem solution
fig, axes = plt.subplots(2, 2, figsize=(14, 10))

# Parameters
a = 1.0  # Cylinder radius
k = 0.1  # Thermal diffusivity
n_terms = 10  # Number of series terms

# Zeros of J_0
zeros = jn_zeros(0, n_terms)

# Initial temperature distribution
def f(r):
    return 1 - r**2  # Parabolic distribution

# Compute Fourier-Bessel coefficients
def compute_coefficients(f, a, zeros, n_terms):
    coeffs = []
    for n in range(n_terms):
        alpha_n = zeros[n]
        
        # Numerator: ∫_0^a r f(r) J_0(α_n r/a) dr
        def integrand_num(r):
            return r * f(r) * jv(0, alpha_n * r / a)
        numerator, _ = quad(integrand_num, 0, a)
        
        # Denominator: ∫_0^a r [J_0(α_n r/a)]^2 dr = (a^2/2) [J_1(α_n)]^2
        denominator = (a**2 / 2) * jv(1, alpha_n)**2
        
        coeffs.append(2 * numerator / (a**2 * jv(1, alpha_n)**2))
    
    return coeffs

coeffs = compute_coefficients(f, a, zeros, n_terms)

# Temperature distribution
def u(r, t, coeffs, zeros, a, k):
    result = np.zeros_like(r)
    for n, (A_n, alpha_n) in enumerate(zip(coeffs, zeros)):
        result += A_n * jv(0, alpha_n * r / a) * np.exp(-k * (alpha_n/a)**2 * t)
    return result

r = np.linspace(0, a, 100)

# Temperature distribution at different times
ax1 = axes[0, 0]
times = [0, 1, 5, 10, 20, 50]
for t in times:
    temp = u(r, t, coeffs, zeros, a, k)
    ax1.plot(r, temp, linewidth=2, label=f't = {t}')

ax1.set_xlabel('Radius r', fontsize=12)
ax1.set_ylabel('Temperature u(r, t)', fontsize=12)
ax1.set_title('Temperature Distribution in Cylinder Over Time', fontsize=14, fontweight='bold')
ax1.legend(fontsize=10)
ax1.grid(True, alpha=0.3)

# Center temperature over time
ax2 = axes[0, 1]
t_array = np.linspace(0, 50, 500)
center_temp = [u(np.array([0.0]), t, coeffs, zeros, a, k)[0] for t in t_array]

ax2.plot(t_array, center_temp, 'b-', linewidth=2)
ax2.set_xlabel('Time t', fontsize=12)
ax2.set_ylabel('Center Temperature u(0, t)', fontsize=12)
ax2.set_title('Temperature at Center vs Time', fontsize=14, fontweight='bold')
ax2.grid(True, alpha=0.3)

# Log scale to see decay
ax3 = axes[1, 0]
ax3.semilogy(t_array, center_temp, 'b-', linewidth=2)
# Theoretically dominated by the first term
theoretical_decay = coeffs[0] * np.exp(-k * (zeros[0]/a)**2 * t_array)
ax3.semilogy(t_array, theoretical_decay, 'r--', linewidth=2, label='First term only')
ax3.set_xlabel('Time t', fontsize=12)
ax3.set_ylabel('Center Temperature (log scale)', fontsize=12)
ax3.set_title('Exponential Decay of Temperature', fontsize=14, fontweight='bold')
ax3.legend(fontsize=10)
ax3.grid(True, alpha=0.3)

# 3D heat map
ax4 = axes[1, 1]
R, T = np.meshgrid(r, t_array[:100])
U = np.zeros_like(R)
for i, t in enumerate(t_array[:100]):
    U[i, :] = u(r, t, coeffs, zeros, a, k)

contour = ax4.contourf(R, T, U, levels=20, cmap='hot')
ax4.set_xlabel('Radius r', fontsize=12)
ax4.set_ylabel('Time t', fontsize=12)
ax4.set_title('Heat Diffusion in Cylinder (Contour)', fontsize=14, fontweight='bold')
plt.colorbar(contour, ax=ax4, label='Temperature')

plt.tight_layout()
plt.savefig('heat_conduction.png', dpi=150, bbox_inches='tight')
plt.show()

print("Fourier-Bessel coefficients (first 5 terms):")
for n, c in enumerate(coeffs[:5]):
    print(f"  A_{n+1} = {c:.6f}")

Generating Functions and Asymptotic Expansions

Generating Function for Bessel FunctionsThis formula connects Bessel functions with complex analysis and is very useful in FM signal analysis.

Generating Function for Legendre Polynomials

Physical significance: This is exactly the expansion of potential energy between two point charges!

If there's a point charge at the origin and another at, the distance between them is.

When:This is the foundation of multipole expansion!

Asymptotic Expansion

As, the behavior of Bessel functions:

Physical significance: Far from the source, cylindrical waves behave as decaying sinusoidal waves.

Connections Between Special Functions

Unified Framework: Sturm-Liouville Theory

All these special functions are eigenfunctions of Sturm-Liouville problems: | Special Function | Equation | Weight| Interval | |-----------------|----------|---------------|----------| ||| 1 || ||| 1 || ||||| ||||| |||||

Common Features

1. Orthogonality: Eigenfunctions for different eigenvalues are orthogonal
2. Completeness: Any "reasonable" function can be expanded
3. Recurrence relations: Simple relations between adjacent orders
4. Generating functions: Compact generating formulas
5. Asymptotic behavior: Predictable behavior at boundaries

Numerical Methods and Implementation

Numerical Computation of Power Series

Direct computation of power series may have numerical issues: 1. Inefficient when many terms are needed 2. Alternating series may have cancellation errors

Better methods: - Use recurrence relations - Use asymptotic formulas (for large arguments) - Tables and interpolation

Special Functions in Python

from scipy import special

# Bessel functions
x = 5.0
print(f"J_0({x}) = {special.jv(0, x)}")
print(f"J_1({x}) = {special.jv(1, x)}")
print(f"Y_0({x}) = {special.yv(0, x)}")

# Zeros of Bessel functions
zeros = special.jn_zeros(0, 5)
print(f"First 5 zeros of J_0: {zeros}")

# Legendre polynomials
print(f"P_3(0.5) = {special.eval_legendre(3, 0.5)}")

# Spherical harmonics
l, m = 2, 1
theta, phi = np.pi/4, np.pi/3
Y = special.sph_harm(m, l, phi, theta)
print(f"Y_{l}^{m}(θ, φ) = {Y}")

# Hermite polynomials
print(f"H_4(1.0) = {special.eval_hermite(4, 1.0)}")

# Airy functions
Ai, Aip, Bi, Bip = special.airy(2.0)
print(f"Ai(2) = {Ai}, Ai'(2) = {Aip}")

Implementing Bessel Functions Yourself

def bessel_j_series(nu, x, n_terms=50):
    """
    Compute Bessel function of the first kind J_ν(x) using series
    
    J_ν(x) = Σ_{m=0}^∞ (-1)^m / (m! Γ(m+ν+1)) * (x/2)^(2m+ν)
    """
    from scipy.special import gamma
    
    result = 0.0
    for m in range(n_terms):
        term = ((-1)**m / (np.math.factorial(m) * gamma(m + nu + 1)) 
                * (x/2)**(2*m + nu))
        result += term
        if abs(term) < 1e-15:  # Convergence check
            break
    return result

# Test
x = 3.0
for nu in [0, 1, 2, 0.5]:
    my_result = bessel_j_series(nu, x)
    scipy_result = special.jv(nu, x)
    print(f"J_{nu}({x}): my implementation = {my_result:.10f}, "
          f"scipy = {scipy_result:.10f}, "
          f"error = {abs(my_result - scipy_result):.2e}")

Hypergeometric Functions: Unifying Special Functions

Gauss Hypergeometric Functionwhereis the Pochhammer symbol.

Amazing fact: Many special functions are special cases of hypergeometric functions!

Confluent Hypergeometric Function

Relation to Hermite polynomials:

Summary: The Power of Series Solutions

Key Points of This Chapter

1. Power series solutions: Expand the solution as, find coefficients via recurrence relations

2. Ordinary points and singular points:

   • Ordinary point: coefficients are analytic, power series solution exists
   • Regular singular point: Frobenius method, generalized power series solution

3. Indicial equation: Algebraic equation determining the power

4. Important special functions:

   • Bessel functions: cylindrical problems
   • Legendre polynomials: spherical problems
   • Hermite polynomials: quantum harmonic oscillator
   • Airy functions: turning point problems

5. Orthogonality and completeness: Expanding arbitrary functions using special functions

6. Unified framework: Sturm-Liouville theory and hypergeometric functions

Preview of Next Chapter

In "Linear Systems of Differential Equations," we will: - Learn about matrix exponentials and fundamental matrices - Analyze the behavior of coupled systems - Study trajectories in phase space - Apply to ecology, circuits, and control theory

Exercises

Basic Problems

1. Use the power series method to solve,,, and verify the answer is.

2. Verify thatis a regular singular point of the equation, and find the indicial equation.

3. Prove thatsatisfies the recurrence relation.

4. Use Rodrigues' formula to computeand.

5. Verify that Hermite polynomials satisfy.

Advanced Problems

6. For the Airy equation:
   • Write the recurrence relation for - Compute the first 8 non-zero coefficients
   • Estimate the radius of convergence
7. Prove that Legendre polynomials satisfy the recurrence relation:$$

(n+1)P_{n+1}(x) = (2n+1)xP_n(x) - nP_{n-1}(x)$J_(x)$satisfies the Bessel differential equation: - Substitute the series into the Bessel equation - Verify the recurrence relation gives correct coefficients

9. Expandonusing Fourier-Bessel series, satisfying. Compute the first 3 coefficients.

Programming Problems

10. Write a function to compute the series expansion of Bessel function, and compare accuracy with scipy.special.jv.

11. Visualize the first 10 energy level wavefunctions of the quantum harmonic oscillator, marking node positions.

12. Simulate circular drumhead vibration:

    • UseandBessel functions
    • Create an animation showing different vibration modes
    • Find the first 5 resonant frequencies

13. Implement multipole expansion:

    • Given a charge distribution, compute its electric potential
    • Expand using Legendre polynomials up to - Compare exact solution with series approximation

Thought Questions

14. Why do so many special functions appear in physics? What is their relationship with the choice of coordinate systems?

15. If the parameterin the Bessel equation is complex, how does the solution change?

16. How does Sturm-Liouville theory unify these seemingly different special functions?

References

1. Arfken, G. B., Weber, H. J., & Harris, F. E. (2012). Mathematical Methods for Physicists. Academic Press.
2. Bender, C. M., & Orszag, S. A. (1999). Advanced Mathematical Methods for Scientists and Engineers. Springer.
3. Abramowitz, M., & Stegun, I. A. (1964). Handbook of Mathematical Functions. Dover.
4. NIST Digital Library of Mathematical Functions: https://dlmf.nist.gov/
5. Watson, G. N. (1995). A Treatise on the Theory of Bessel Functions. Cambridge University Press.

---

Next Chapter: Linear Systems of Differential Equations →

---

This is Chapter 5 of the "World of Ordinary Differential Equations" series.