PDEs & Separation of Variables

Partial differential equations are the language of continuum physics — from heat conduction to wave propagation to quantum mechanics. The method of separation of variables reduces PDEs to systems of ODEs by exploiting coordinate symmetry, producing eigenfunction expansions that form the backbone of mathematical physics.

1. Classification of Second-Order PDEs

The most general second-order linear PDE in two independent variables is:

$$A\frac{\partial^2 u}{\partial x^2} + 2B\frac{\partial^2 u}{\partial x\,\partial y} + C\frac{\partial^2 u}{\partial y^2} + D\frac{\partial u}{\partial x} + E\frac{\partial u}{\partial y} + Fu = G$$

The classification depends on the discriminant $\Delta = B^2 - AC$, analogous to the classification of conic sections:

  • Elliptic ($\Delta < 0$): Laplace/Poisson equation. Boundary value problems. Steady-state phenomena.
  • Parabolic ($\Delta = 0$): Heat/diffusion equation. Initial-boundary value problems. Irreversible processes.
  • Hyperbolic ($\Delta > 0$): Wave equation. Initial value problems. Wave propagation with finite speed.

1.1 Derivation of Classification via Characteristics

The classification arises from the theory of characteristics. Consider the principal part $Au_{xx} + 2Bu_{xy} + Cu_{yy} = 0$. Along a curve $y = y(x)$, we seek curves where the second derivatives are not uniquely determined by the equation and data along the curve. Setting $\xi = dy/dx$, the characteristic equation is:

$$A\xi^2 - 2B\xi + C = 0 \implies \xi = \frac{B \pm \sqrt{B^2 - AC}}{A}$$

The nature of the roots determines the type:

  • $B^2 - AC > 0$: Two real characteristics (hyperbolic) — information propagates along these curves
  • $B^2 - AC = 0$: One repeated real characteristic (parabolic) — one-sided propagation
  • $B^2 - AC < 0$: No real characteristics (elliptic) — information spreads everywhere simultaneously

1.2 The Canonical Equations of Physics

Laplace's equation (elliptic): $A = C = 1$, $B = 0$, so $\Delta = -1 < 0$:

$$\nabla^2 u = \frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2} = 0$$

Heat equation (parabolic): With $y \to t$, $A = \alpha$, $B = 0$, $C = 0$, $E = -1$:

$$\frac{\partial u}{\partial t} = \alpha \frac{\partial^2 u}{\partial x^2}$$

Wave equation (hyperbolic): $A = c^2$, $B = 0$, $C = -1$ (with $y \to t$), so $\Delta = c^2 > 0$:

$$\frac{\partial^2 u}{\partial t^2} = c^2 \frac{\partial^2 u}{\partial x^2}$$

2. Separation of Variables: General Method

The method of separation of variables assumes a product solution and reduces the PDE to ODEs. The key steps are:

  1. Assume $u(\mathbf{r}, t) = F_1(x_1) F_2(x_2) \cdots$ as a product of functions of single variables.
  2. Substitute into the PDE and divide to isolate terms depending on different variables.
  3. Since each separated term depends on a different variable, each must equal a constant (the separation constant).
  4. Solve the resulting ODEs with appropriate boundary conditions to find eigenvalues and eigenfunctions.
  5. Form the general solution as a superposition (sum or integral) of product solutions.
  6. Determine coefficients from initial conditions using orthogonality.

2.1 Why Separation Works: Symmetry and Coordinate Systems

Separation of variables succeeds when the PDE, domain, and boundary conditions all respect the symmetry of a coordinate system. The Laplacian $\nabla^2$ is separable in exactly 11 orthogonal coordinate systems in 3D (Cartesian, cylindrical, spherical, elliptic cylindrical, parabolic cylindrical, prolate/oblate spheroidal, paraboloidal, ellipsoidal, conical, and bispherical). We focus on the three most important.

3. Separation in Cartesian Coordinates

3.1 The 1D Heat Equation

Consider heat conduction in a rod of length $L$ with endpoints held at zero temperature:

$$\frac{\partial u}{\partial t} = \alpha \frac{\partial^2 u}{\partial x^2}, \quad u(0,t) = u(L,t) = 0, \quad u(x,0) = f(x)$$

Step 1: Assume $u(x,t) = X(x)T(t)$. Substituting:

$$X(x)T'(t) = \alpha X''(x)T(t) \implies \frac{T'(t)}{\alpha T(t)} = \frac{X''(x)}{X(x)} = -\lambda$$

Step 2: The spatial ODE with boundary conditions:

$$X'' + \lambda X = 0, \quad X(0) = 0, \quad X(L) = 0$$

This is a Sturm-Liouville eigenvalue problem. For $\lambda > 0$, the general solution is$X = A\sin(\sqrt{\lambda}\,x) + B\cos(\sqrt{\lambda}\,x)$. The boundary condition$X(0) = 0$ gives $B = 0$, and $X(L) = 0$ requires$\sin(\sqrt{\lambda}\,L) = 0$, giving eigenvalues:

$$\lambda_n = \left(\frac{n\pi}{L}\right)^2, \quad X_n(x) = \sin\!\left(\frac{n\pi x}{L}\right), \quad n = 1, 2, 3, \ldots$$

Step 3: The temporal ODE $T' + \alpha\lambda_n T = 0$ gives$T_n(t) = e^{-\alpha\lambda_n t}$. The general solution is:

$$u(x,t) = \sum_{n=1}^{\infty} b_n \sin\!\left(\frac{n\pi x}{L}\right) e^{-\alpha(n\pi/L)^2 t}$$

where the coefficients are determined by the initial condition using orthogonality:

$$b_n = \frac{2}{L}\int_0^L f(x)\sin\!\left(\frac{n\pi x}{L}\right)dx$$

3.2 The 1D Wave Equation

For a vibrating string of length $L$ fixed at both ends:

$$\frac{\partial^2 u}{\partial t^2} = c^2 \frac{\partial^2 u}{\partial x^2}, \quad u(0,t) = u(L,t) = 0$$

Separating $u = X(x)T(t)$ gives the same spatial eigenvalue problem, and the temporal equation becomes $T'' + c^2\lambda_n T = 0$ with solutions $T_n = A_n\cos(\omega_n t) + B_n\sin(\omega_n t)$where $\omega_n = cn\pi/L$:

$$u(x,t) = \sum_{n=1}^{\infty} \sin\!\left(\frac{n\pi x}{L}\right)\!\left[A_n\cos\!\left(\frac{n\pi ct}{L}\right) + B_n\sin\!\left(\frac{n\pi ct}{L}\right)\right]$$

The coefficients $A_n$ and $B_n$ are fixed by the initial displacement$u(x,0)$ and initial velocity $u_t(x,0)$. These are the normal modesof the string, each oscillating at a frequency that is an integer multiple of the fundamental$\omega_1 = c\pi/L$.

3.3 Laplace's Equation in a Rectangle

Consider $\nabla^2 u = 0$ in $0 < x < a$, $0 < y < b$ with$u = 0$ on three sides and $u(x, b) = f(x)$. Separating $u = X(x)Y(y)$:

$$\frac{X''}{X} = -\frac{Y''}{Y} = -\lambda$$

The $x$-equation with homogeneous BCs gives $X_n = \sin(n\pi x/a)$,$\lambda_n = (n\pi/a)^2$. The $y$-equation $Y'' - \lambda_n Y = 0$with $Y(0) = 0$ gives $Y_n = \sinh(n\pi y/a)$:

$$u(x,y) = \sum_{n=1}^{\infty} c_n \sin\!\left(\frac{n\pi x}{a}\right) \sinh\!\left(\frac{n\pi y}{a}\right)$$

4. Separation in Cylindrical Coordinates

In cylindrical coordinates $(r, \phi, z)$, the Laplacian is:

$$\nabla^2 u = \frac{1}{r}\frac{\partial}{\partial r}\!\left(r\frac{\partial u}{\partial r}\right) + \frac{1}{r^2}\frac{\partial^2 u}{\partial \phi^2} + \frac{\partial^2 u}{\partial z^2}$$

Setting $u = R(r)\Phi(\phi)Z(z)$ and separating gives three ODEs. The azimuthal equation$\Phi'' + m^2\Phi = 0$ with periodicity $\Phi(\phi + 2\pi) = \Phi(\phi)$gives $m = 0, 1, 2, \ldots$ and $\Phi_m = e^{\pm im\phi}$. The$z$-equation gives exponentials or trigonometric functions depending on the separation constant $k$.

4.1 Bessel's Equation

The radial equation becomes Bessel's equation:

$$r^2 R'' + rR' + (k^2 r^2 - m^2)R = 0$$

With the substitution $\rho = kr$, this is Bessel's equation of order $m$. The solutions are Bessel functions of the first kind $J_m(kr)$ (regular at $r=0$) and second kind $Y_m(kr)$ (singular at $r=0$).

For a cylinder of radius $a$ with $u(a,\phi,z) = 0$, the boundary condition$J_m(ka) = 0$ gives discrete eigenvalues $k_{mn} = j_{mn}/a$ where$j_{mn}$ is the $n$-th zero of $J_m$.

4.2 Heat Equation in a Cylinder

For the heat equation in a cylinder of radius $a$ with no $z$-dependence and azimuthal symmetry ($m = 0$):

$$u(r,t) = \sum_{n=1}^{\infty} c_n J_0\!\left(\frac{j_{0n}\,r}{a}\right) e^{-\alpha(j_{0n}/a)^2 t}$$

where $j_{0n}$ are the zeros of $J_0$ and the coefficients use the orthogonality of Bessel functions with weight $r$:

$$c_n = \frac{2}{a^2 [J_1(j_{0n})]^2} \int_0^a r\, f(r)\, J_0\!\left(\frac{j_{0n}\,r}{a}\right) dr$$

5. Separation in Spherical Coordinates

In spherical coordinates $(r, \theta, \phi)$, the Laplacian is:

$$\nabla^2 u = \frac{1}{r^2}\frac{\partial}{\partial r}\!\left(r^2\frac{\partial u}{\partial r}\right) + \frac{1}{r^2\sin\theta}\frac{\partial}{\partial \theta}\!\left(\sin\theta\frac{\partial u}{\partial \theta}\right) + \frac{1}{r^2\sin^2\theta}\frac{\partial^2 u}{\partial \phi^2}$$

Setting $u = R(r)\Theta(\theta)\Phi(\phi)$ and separating:

5.1 The Angular Equations

The azimuthal equation gives $\Phi_m(\phi) = e^{im\phi}$ with $m \in \mathbb{Z}$. The polar equation, with the substitution $x = \cos\theta$, becomes the associated Legendre equation:

$$\frac{d}{dx}\!\left[(1-x^2)\frac{dP}{dx}\right] + \left[\ell(\ell+1) - \frac{m^2}{1-x^2}\right]P = 0$$

Regularity at $x = \pm 1$ (the poles) requires $\ell = 0, 1, 2, \ldots$and $|m| \le \ell$. The solutions are the associated Legendre functions$P_\ell^m(\cos\theta)$. Combined with the azimuthal part, we get the spherical harmonics:

$$Y_\ell^m(\theta, \phi) = \sqrt{\frac{(2\ell+1)}{4\pi}\frac{(\ell-m)!}{(\ell+m)!}}\; P_\ell^m(\cos\theta)\, e^{im\phi}$$

These form a complete orthonormal set on the sphere:

$$\int_0^{2\pi}\!\int_0^{\pi} Y_{\ell'}^{m'*}(\theta,\phi)\, Y_\ell^m(\theta,\phi)\,\sin\theta\,d\theta\,d\phi = \delta_{\ell\ell'}\delta_{mm'}$$

5.2 The Radial Equation

The radial equation for Laplace's equation is:

$$\frac{d}{dr}\!\left(r^2\frac{dR}{dr}\right) = \ell(\ell+1)R \implies R(r) = A_\ell r^\ell + B_\ell r^{-(\ell+1)}$$

The general solution to Laplace's equation in spherical coordinates is therefore:

$$u(r,\theta,\phi) = \sum_{\ell=0}^{\infty}\sum_{m=-\ell}^{\ell} \left(A_{\ell m} r^\ell + B_{\ell m} r^{-(\ell+1)}\right) Y_\ell^m(\theta,\phi)$$

5.3 Example: Potential of a Sphere

For the potential outside a sphere of radius $a$ with prescribed surface potential$u(a,\theta,\phi) = f(\theta,\phi)$, regularity at infinity requires $A_{\ell m} = 0$:

$$u(r,\theta,\phi) = \sum_{\ell=0}^{\infty}\sum_{m=-\ell}^{\ell} f_{\ell m} \left(\frac{a}{r}\right)^{\ell+1} Y_\ell^m(\theta,\phi)$$

where $f_{\ell m} = \int Y_\ell^{m*} f\, d\Omega$. This is the multipole expansion fundamental to electrostatics, gravitational theory, and quantum mechanics.

6. The Helmholtz Equation and Normal Modes

The time-independent part of both the wave equation and the Schrodinger equation leads to the Helmholtz equation:

$$\nabla^2 u + k^2 u = 0$$

This arises whenever we separate time from a wave-like equation. For example, if$u(\mathbf{r}, t) = \psi(\mathbf{r})e^{-i\omega t}$ in the wave equation, then$\nabla^2\psi + ({\omega^2}/{c^2})\psi = 0$.

In quantum mechanics, the time-independent Schrodinger equation$-\frac{\hbar^2}{2m}\nabla^2\psi + V\psi = E\psi$ is a Helmholtz-type equation where $k^2 = 2m(E-V)/\hbar^2$.

Normal modes of a drum: The Helmholtz equation on a circular membrane of radius$a$ with $\psi(a,\phi) = 0$ gives modes:

$$\psi_{mn}(r,\phi) = J_m\!\left(\frac{j_{mn}\,r}{a}\right)\begin{cases}\cos(m\phi) \\ \sin(m\phi)\end{cases}, \quad k_{mn} = \frac{j_{mn}}{a}$$

The frequencies $\omega_{mn} = c\, j_{mn}/a$ are not integer multiples of the fundamental (unlike a string), which is why drums produce less harmonious tones than strings.

7. Well-Posedness and Uniqueness Theorems

A problem is well-posed (in the sense of Hadamard) if: (1) a solution exists, (2) the solution is unique, and (3) the solution depends continuously on the data.

Maximum principle for harmonic functions: If $u$ is harmonic in a bounded domain $D$ and continuous on $\bar{D}$, then $u$ attains its maximum and minimum on the boundary $\partial D$. This immediately implies uniqueness for the Dirichlet problem (if two solutions agree on the boundary, their difference is harmonic with zero boundary data, hence zero everywhere).

Energy methods for the wave equation: Define the energy:

$$E(t) = \frac{1}{2}\int_\Omega\!\left[\left(\frac{\partial u}{\partial t}\right)^{\!2} + c^2|\nabla u|^2\right]d\mathbf{r}$$

For homogeneous boundary conditions, $dE/dt = 0$ (energy conservation). If two solutions share the same initial and boundary data, their difference has zero energy, hence is identically zero.

8. d'Alembert's Solution and Method of Characteristics

For the wave equation on an infinite domain, there is an elegant closed-form solution that reveals the wave nature directly. Introducing characteristic coordinates$\xi = x - ct$ and $\eta = x + ct$:

$$\frac{\partial^2 u}{\partial t^2} = c^2\frac{\partial^2 u}{\partial x^2} \implies \frac{\partial^2 u}{\partial \xi\,\partial \eta} = 0$$

The general solution is $u = F(\xi) + G(\eta) = F(x - ct) + G(x + ct)$, a superposition of right-traveling and left-traveling waves. Imposing initial conditions$u(x,0) = f(x)$ and $u_t(x,0) = g(x)$:

$$u(x,t) = \frac{1}{2}\left[f(x-ct) + f(x+ct)\right] + \frac{1}{2c}\int_{x-ct}^{x+ct} g(s)\,ds$$

This is d'Alembert's formula. It shows that the solution at any point$(x, t)$ depends only on initial data in the interval $[x - ct, x + ct]$ — the domain of dependence. Information propagates at finite speed $c$, in stark contrast to the heat equation where disturbances propagate instantaneously.

Physical interpretation: The initial displacement $f(x)$ splits into two halves, each propagating in opposite directions. The initial velocity $g(x)$generates waves that spread outward from each point. Characteristics$x \pm ct = \text{const}$ are the paths along which wave information travels.

For the heat equation, the characteristic analysis yields only one family of characteristics ($t = \text{const}$), reflecting the fact that the equation is parabolic and information propagates in only one direction in time.

9. Boundary Conditions and Their Physical Meaning

The choice of boundary conditions profoundly affects the solution and its physical interpretation:

9.1 Dirichlet Conditions

$u = f$ on the boundary. Physically: fixed temperature (heat), fixed displacement (vibrations), specified potential (electrostatics). The eigenvalues for $-\nabla^2$ on a domain with Dirichlet conditions are all positive.

9.2 Neumann Conditions

$\partial u/\partial n = g$ on the boundary (normal derivative). Physically: insulated boundary (heat), free endpoint (vibrations), specified charge density (electrostatics). Neumann conditions allow a zero eigenvalue (constant eigenfunction), and a solution to the Poisson equation exists only if the source satisfies a compatibility condition:

$$\int_\Omega f\,d\mathbf{r} = \int_{\partial\Omega} g\,dS$$

9.3 Robin (Mixed) Conditions

$\alpha u + \beta \partial u/\partial n = g$. Physically: convective heat transfer (Newton's law of cooling), impedance boundary in acoustics, partially absorbing boundaries in quantum mechanics.

9.4 Periodic Conditions

$u(a) = u(b)$, $u'(a) = u'(b)$. These arise naturally in problems with angular coordinates (the azimuthal angle $\phi$ in cylindrical or spherical coordinates). Periodic boundary conditions can admit degenerate eigenvalues (e.g., $\cos(n\phi)$ and$\sin(n\phi)$ share the same eigenvalue).

8. Python Simulation: Separation of Variables

This simulation numerically solves the heat equation and wave equation using separation of variables, verifies eigenfunction orthogonality, and demonstrates convergence of the Fourier series solutions.

Python
script.py241 lines

Click Run to execute the Python code

Code will be executed with Python 3 on the server