ODE Chapter 5: Power Series and Special Functions

Some ODEs have no solutions in terms of familiar functions. The Bessel equation, the Legendre equation, the Airy equation – all arise naturally in physics (heat conduction in cylinders, gravitational fields of planets, quantum tunneling). Their solutions define entirely new functions. This chapter shows you how to find them using power series, why the Frobenius extension is forced upon us at singular points, and why the same handful of “special functions” keeps appearing across physics and engineering.

What You Will Learn

• Power series solutions at ordinary points and the role of the radius of convergence
• Regular singular points and the Frobenius method
• Bessel functions: vibrating drums and cylindrical heat conduction
• Legendre polynomials: spherical harmonics and multipole expansions
• Hermite polynomials: the quantum harmonic oscillator
• Airy functions: turning points in quantum mechanics
• The Sturm-Liouville framework that unifies all of the above

Prerequisites

• Taylor series and radius of convergence
• Chapters 1-3: second-order linear ODEs
• Basic complex numbers (we will use the complex plane to read off the radius of convergence)

1. Why Power Series?

This is the Bessel equation. None of our previous methods (characteristic equations for constant coefficients, integrating factors, variation of parameters) reach it because the coefficients are functions of $x$, not constants. Worse, the equation is singular at $x = 0$ – precisely where the cylinder’s axis lies and precisely where the answer must be physically well-defined.

The strategy. Postulate a series ansatz, plug it in, and let the differential equation generate its own coefficients through a recurrence. We will need two flavours of series:

1. Plain Taylor series $\sum a_k (x - x_0)^k$ at ordinary points.
2. Frobenius series $(x - x_0)^r \sum a_k (x - x_0)^k$ at regular singular points, where the exponent $r$ is itself determined by the equation.

The same machinery unlocks Bessel, Legendre, Hermite, Airy, Chebyshev, Laguerre, hypergeometric – the entire pantheon of classical special functions.

2. Power Series at Ordinary Points

2.1 The method

A point $x_0$ is ordinary if both $P$ and $Q$ are analytic there. The recipe is then:

1. Verify $x_0$ is ordinary.
2. Postulate $y = \sum_{k=0}^{\infty} a_k (x - x_0)^k$.
3. Differentiate term by term and substitute into the ODE.
4. Collect coefficients of each power $(x - x_0)^k$ and equate to zero.
5. Solve the resulting recurrence for $a_k$ in terms of $a_0$ and $a_1$ (the two free parameters that match the two initial conditions).

A theorem of Fuchs guarantees that the resulting series converges in any open disk centred at $x_0$ that contains no singularity of $P$ or $Q$ – and the radius extends all the way out to the nearest singularity in the complex plane. This is why complex analysis matters even when the ODE itself is real-valued.

2.2 Warm-up: rediscovering $e^x$

The series method “rediscovered” the exponential. The lesson: even when we already know the answer, the recurrence is a mechanical procedure that cannot fail.

2.3 The Airy equation

This appears whenever a wave (light, quantum probability, water) meets a turning point where the local wavenumber smoothly changes from real to imaginary – the prototype example is quantum tunnelling near a linear potential barrier.

The coefficients split into two independent chains starting from $a_0$ and $a_1$, and the resulting two linearly independent solutions are the Airy functions $\text{Ai}(x)$ and $\text{Bi}(x)$. Airy decays exponentially for $x > 0$ (the classically forbidden region) and oscillates for $x < 0$ (the classically allowed region) – a crisp visualisation of quantum tunnelling.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16

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

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

plt.figure(figsize=(10, 5))
plt.plot(x, Ai, 'b-', linewidth=2, label='Ai(x)')
plt.plot(x, Bi, 'r-', linewidth=2, label='Bi(x)')
plt.axhline(y=0, color='k', linewidth=0.5)
plt.xlabel('x'); plt.ylabel('y')
plt.title("Airy Functions: Solutions of y'' - xy = 0")
plt.legend(); plt.ylim(-0.6, 1.2)
plt.grid(True, alpha=0.3); plt.tight_layout()
plt.show()

3. Regular Singular Points and the Frobenius Method

3.1 Why a plain Taylor series can fail

By inspection, $y = \sqrt{x}$ is a solution – but $\sqrt{x}$ has a vertical tangent at the origin and is not analytic there, so no series of the form $\sum a_k x^k$ can ever reproduce it. Some “fractional power” of $x$ has to enter the ansatz.

3.2 Regular singular points

are both analytic at $x_0$. Roughly: the singularity is at most a simple pole of $P$ and at most a double pole of $Q$.

3.3 The Frobenius ansatz

where $p_0 = \lim_{x \to x_0}(x-x_0) P(x)$ and $q_0 = \lim_{x \to x_0}(x-x_0)^2 Q(x)$. Its two roots $r_1 \geq r_2$ govern how to assemble two independent solutions:

The recurrence on $a_k$, obtained from the higher-order coefficient matches, is then unwound term by term.

4. Bessel Functions

4.1 The Bessel equation

The three places this equation appears in physics:

• Vibrating circular drumheads – the radial part of $\nabla^2 u + k^2 u = 0$ in polar coordinates.
• Cylindrical waveguides – electromagnetic modes inside a metal pipe.
• Heat diffusion in a cylinder – transient cooling of a metal rod.

The point $x = 0$ is a regular singular point with $p_0 = 1$ and $q_0 = -n^2$, so the indicial equation $r^2 - n^2 = 0$ gives $r = \pm n$.

4.2 The Bessel function of the first kind

Three properties we will use repeatedly:

• Recurrence: $\displaystyle J_{n-1}(x) + J_{n+1}(x) = \frac{2n}{x}\,J_n(x)$ – lets us slide between orders.
• Asymptotic behaviour: $\displaystyle J_n(x) \approx \sqrt{\frac{2}{\pi x}}\cos\!\left(x - \frac{n\pi}{2} - \frac{\pi}{4}\right)$ as $x \to \infty$ – so $J_n$ behaves eventually like a damped cosine.
• Zeros: the positive zeros $j_{n,1} < j_{n,2} < \dots$ of $J_n$ are spaced approximately $\pi$ apart and quantise the resonant frequencies of a circular membrane.

4.3 From Bessel zeros to drumhead modes

where $j_{n,m}$ is the $m$-th positive zero of $J_n$. The integer $n$ counts the number of angular nodal lines, and $m$ counts the number of radial nodal circles (including the boundary).

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21

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

x = np.linspace(0.01, 20, 1000)
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 4))

for n in range(5):
    ax1.plot(x, jv(n, x), linewidth=2, label=f'$J_{n}(x)$')
ax1.axhline(0, color='k', linewidth=0.5)
ax1.set_xlabel('x'); ax1.set_title('Bessel Functions of the First Kind')
ax1.legend(fontsize=9); ax1.set_ylim(-0.5, 1.1); ax1.grid(True, alpha=0.3)

for n in range(4):
    zeros = jn_zeros(n, 5)
    ax2.scatter(zeros, [n]*5, s=80, label=f'$J_{n}$ zeros')
ax2.set_xlabel('x'); ax2.set_ylabel('Order n')
ax2.set_title('Zeros of Bessel Functions (Drum Modes)')
ax2.legend(fontsize=9); ax2.grid(True, alpha=0.3)

plt.tight_layout(); plt.show()

4.4 Bessel function of the second kind

For the negative root $r = -n$ when $n$ is a non-negative integer, the recurrence breaks down: a logarithmic term is forced. The resulting linearly independent partner is the Bessel function of the second kind $Y_n(x)$, which behaves as $Y_n(x) \sim -(2/\pi)(n-1)!\,(x/2)^{-n}$ near the origin and so blows up there. Physical solutions on a domain that includes the axis (a solid drum, a solid wire) must therefore use only $J_n$; problems on an annular domain (a coaxial cable) need both $J_n$ and $Y_n$.

5. Legendre Polynomials

5.1 The Legendre equation

This appears every time you separate variables in spherical coordinates: the angular part of the Laplacian, hydrogen-atom wavefunctions, gravitational and electrostatic multipole expansions, antenna radiation patterns. Here $x$ stands in for $\cos\theta$, so the natural domain is $[-1,1]$.

The points $x = \pm 1$ are regular singular points (where the leading coefficient $1-x^2$ vanishes), while $x = 0$ is ordinary.

5.2 Why polynomials?

When $n$ is a non-negative integer, the numerator vanishes at $k = n$ and the series truncates to a polynomial. The other linearly independent solution remains an infinite series and diverges at $x = \pm 1$, so demanding regularity at the poles of the sphere selects the polynomial solutions and quantises the parameter $n$.

Normalised so that $P_n(1) = 1$:

5.3 Rodrigues’ formula and orthogonality

Any reasonable function on $[-1,1]$ admits a Legendre expansion $f(x) = \sum_{n=0}^{\infty} c_n P_n(x)$ with $c_n = \frac{2n+1}{2}\int_{-1}^{1} f(x) P_n(x)\,dx$ – the basis of the multipole expansion in electromagnetism and gravitation.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13

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

x = np.linspace(-1, 1, 500)
plt.figure(figsize=(8, 5))
for n in range(6):
    Pn = legendre(n)
    plt.plot(x, Pn(x), linewidth=2, label=f'$P_{n}(x)$')
plt.axhline(0, color='k', linewidth=0.5)
plt.xlabel('x'); plt.title('Legendre Polynomials')
plt.legend(fontsize=9); plt.grid(True, alpha=0.3)
plt.tight_layout(); plt.show()

6. Hermite Polynomials and the Quantum Harmonic Oscillator

6.1 The Hermite equation

They satisfy the orthogonality $\int_{-\infty}^{\infty} H_m(x) H_n(x) e^{-x^2} dx = \sqrt{\pi}\,2^n n!\,\delta_{mn}$ – with the Gaussian weight that will become the wavefunction envelope in a moment.

6.2 The quantum harmonic oscillator

Two physical lessons fall out for free:

• Energy quantisation. Only the integer values $n = 0, 1, 2, \ldots$ give normalisable solutions, so the energy ladder is discrete.
• Zero-point energy. Even the ground state has $E_0 = \tfrac{1}{2}\hbar\omega \neq 0$ – a quantum particle in a parabolic well cannot sit perfectly still.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18

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

x = np.linspace(-4, 4, 500)
plt.figure(figsize=(8, 6))
for n in range(5):
    Hn = hermite(n)
    norm = 1 / np.sqrt(2**n * factorial(n)) * (1/np.pi)**0.25
    psi = norm * Hn(x) * np.exp(-x**2 / 2)
    plt.plot(x, psi + n, linewidth=2, label=f'n={n}')
    plt.fill_between(x, n, psi + n, alpha=0.2)
    plt.axhline(n, color='gray', linewidth=0.5, linestyle='--')

plt.xlabel('Position x'); plt.ylabel('Wavefunction + energy offset')
plt.title('Quantum Harmonic Oscillator Wavefunctions')
plt.grid(True, alpha=0.3); plt.tight_layout(); plt.show()

7. The Sturm-Liouville Unifying Framework

with appropriate boundary conditions. Different choices of $p$, $q$, $w$, and the interval reproduce each special function:

Three properties are inherited by every entry:

1. Real eigenvalues – the operator is symmetric under the weight $w$.
2. Orthogonality – $\int y_m y_n\,w\,dx = 0$ for $m \neq n$.
3. Completeness – generic functions admit a generalised Fourier series in the eigenfunctions.

This is why “expand in spherical harmonics”, “expand in Bessel modes”, “expand in Hermite states” all feel the same: they are the same theorem with different weights.

8. Python: Using Special Functions

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16

from scipy import special
import numpy as np

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

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

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

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

The figures in this article were produced by scripts/figures/ode/05-power-series.py; rerun it after pulling to regenerate every PNG into both the EN and ZH asset folders.

Summary

Exercises

Basic

1. Use power series to solve $y'' + y = 0$ near $x = 0$. Verify you recover $\sin x$ and $\cos x$.
2. Show that $x = 0$ is a regular singular point of $2x^2 y'' + 3x y' - (1 + x) y = 0$ and find its indicial equation.
3. Compute $P_4(x)$ and $P_5(x)$ from Rodrigues’ formula and check the recurrence $(n+1) P_{n+1}(x) = (2n+1) x P_n(x) - n P_{n-1}(x)$.

Advanced

4. For the Airy equation $y'' = xy$, write the recurrence and compute the first 8 non-zero coefficients of each independent solution.
5. Prove the Bessel recurrence $J_{n-1}(x) + J_{n+1}(x) = (2n/x) J_n(x)$ from the series for $J_n$.
6. Expand $f(r) = 1 - r^2$ on $[0, 1]$ as a Fourier-Bessel series in $J_0$ and verify the first few coefficients numerically.

Programming

7. Implement the series for $J_n(x)$ from scratch and compare with scipy.special.jv over $x \in [0, 30]$ for $n = 0, 1, 2$.
8. Re-create the drumhead figure for modes $(n, m) \in \{(0,1),(0,2),(1,1),(2,1),(3,1)\}$ and animate one period in time.
9. Plot the first 10 quantum harmonic oscillator probability densities $|\psi_n|^2$ on top of the parabolic potential, marking the classical turning points $x = \pm\sqrt{2n+1}$.

References

• Arfken, Weber & Harris, Mathematical Methods for Physicists, Academic Press (2012)
• Bender & Orszag, Advanced Mathematical Methods, Springer (1999)
• NIST Digital Library of Mathematical Functions: dlmf.nist.gov ↗

Series Navigation