Separation of Variables Method

Learning Objectives

By the end of this section, you will be able to:

1. Explain the key assumption behind separation of variables: that the solution is a product $u(x,t) = X(x) \cdot T(t)$
2. Execute the method step by step — substitute, separate, set equal to a constant, and solve two ODEs
3. Solve eigenvalue problems arising from boundary conditions and identify the eigenfunctions
4. Apply superposition to construct general solutions from individual modes
5. Compute Fourier coefficients using orthogonality of eigenfunctions
6. Contrast how the method produces exponential decay for the heat equation versus oscillation for the wave equation
7. Connect separation of variables to spectral methods and Fourier Neural Operators in modern machine learning

The Big Picture: From One PDE to Many ODEs

"The art of doing mathematics is finding that special case which contains all the germs of generality." — David Hilbert

In the previous sections, we learned what partial differential equations are, how they are classified, and what boundary and initial conditions mean. Now comes the central question: how do we actually solve a PDE?

Partial differential equations are fundamentally harder than ordinary differential equations. An ODE involves derivatives with respect to one variable; a PDE involves derivatives with respect to multiple variables simultaneously. The solution is not a curve but an entire surface (or higher-dimensional object). We cannot simply "integrate both sides."

The method of separation of variables, developed by Joseph Fourier in the early 1800s while studying heat conduction, offers a profound insight: we can break a complicated PDE into a collection of simpler ODEs. The trick is to assume that the solution has a special multiplicative structure and then show that this assumption leads to a complete family of solutions.

Historical Context

Joseph Fourier's 1822 masterwork Théorie Analytique de la Chaleur (Analytical Theory of Heat) introduced both the heat equation and the method of separation of variables. His key realization was that any reasonable function could be represented as an infinite sum of sines and cosines — what we now call a Fourier series. This was initially controversial (even Lagrange and Laplace were skeptical), but it revolutionized both mathematics and physics.

The method was later applied to the wave equation by Daniel Bernoulli and to Laplace's equation by many others. Today, it remains one of the most important analytical tools in mathematical physics and engineering, and its ideas are the foundation of spectral methods in computational science and Fourier Neural Operators in machine learning.

The Core Idea: Products Become Sums

The entire method rests on a single, inspired guess:

Why would this work? Not every function of two variables is a product of two functions of one variable. The function $u(x,t) = x + t$ cannot be written as $X(x) \cdot T(t)$. But the key insight is that we don't need every solution to have this form — we only need enough of them that we can build the general solution by adding them together (superposition).

The Method Step by Step

Let us develop the method using the one-dimensional heat equation as our prototype. The same steps apply (with minor modifications) to the wave equation, Laplace's equation, and many other PDEs.

Step 1: State the Problem

We want to solve the heat equation on a rod of length $L$ with both ends held at zero temperature:

Step 2: Assume a Product Solution

Substitute $u(x, t) = X(x) \cdot T(t)$ into the PDE:

Here $T'(t)$ means $dT/dt$ and $X''(x)$ means $d^2X/dx^2$. Note how the product structure allows partial derivatives to act on only one factor.

Step 3: Divide to Separate

Divide both sides by $X(x) \cdot T(t)$ (assuming neither is zero):

Now comes the critical observation. The left side depends only on $t$. The right side depends only on $x$. Yet they are equal for all values of both $x$ and $t$. The only way this can happen is if both sides equal the same constant.

Step 4: Introduce the Separation Constant

We call this constant $-\lambda$ (the minus sign and the symbol $\lambda$ — lambda — are chosen by convention for convenience):

We have converted one PDE into two ODEs! Each is a standard ODE that we know how to solve. The parameter $\lambda$ connects them: it is the separation constant (which will turn out to be an eigenvalue).

Step 5: Apply Boundary Conditions

The boundary conditions $u(0, t) = 0$ and $u(L, t) = 0$ translate directly to the spatial part (since $T(t)$ is not identically zero):

This transforms the spatial ODE into an eigenvalue problem: find the values of $\lambda$ for which the ODE $X'' + \lambda X = 0$ has a nonzero solution satisfying both boundary conditions.

Example: The Heat Equation on a Finite Rod

Let us carry out the full solution for the heat equation with Dirichlet boundary conditions (both ends at zero temperature).

Solving the Spatial ODE

The spatial ODE is:

The general solution depends on the sign of $\lambda$. We consider three cases:

Therefore, the eigenvalues and eigenfunctions are:

Solving the Temporal ODE

For each eigenvalue $\lambda_n$, the temporal ODE is:

This is a first-order linear ODE with solution:

Each mode decays exponentially, with higher modes (larger $n$) decaying much faster because the decay rate is proportional to $n^2$.

Individual Mode Solutions

Combining spatial and temporal parts, each mode (or normal mode) of the heat equation is:

Eigenvalue Problems: The Heart of the Method

The spatial ODE together with its boundary conditions forms an eigenvalue problem. This is one of the most important concepts in mathematical physics and has deep connections to linear algebra.

The analogy with linear algebra is powerful. In matrix algebra, we have $Av = \lambda v$ where $A$ is a matrix and $v$ is an eigenvector. Here, the differential operator $-d^2/dx^2$ plays the role of the matrix, and the eigenfunctions $\sin(n\pi x/L)$ play the role of eigenvectors.

Different boundary conditions produce different eigenfunctions:

Superposition and Fourier Series

Each individual mode $u_n(x,t) = \sin(n\pi x/L) \, e^{-\alpha(n\pi/L)^2 t}$ solves the PDE and satisfies the boundary conditions. But there are infinitely many such solutions (one for each $n$). Which one is the actual solution?

The answer uses superposition: because the heat equation is linear, any linear combination of solutions is also a solution. Therefore, the general solution is:

The coefficients $b_n$ are determined by the initial condition $u(x, 0) = f(x)$. At $t = 0$, all the exponentials equal 1, so:

This is a Fourier sine series! The coefficients are computed using orthogonality: multiply both sides by $\sin(m\pi x/L)$ and integrate from 0 to $L$. Because of orthogonality, all terms vanish except $n = m$:

Why Orthogonality Matters

Without orthogonality, finding the coefficients $b_n$ would require solving an infinite system of coupled equations — essentially impossible. Orthogonality allows us to compute each coefficient independently, using a single integral. This is analogous to how orthogonal coordinate axes let us read off each component of a vector independently.

Convergence and the Complete Solution

The full solution to the heat equation is then:

This formula tells a beautiful physical story: the initial temperature distribution $f(x)$ is decomposed into sinusoidal modes, each of which decays independently at a rate proportional to $n^2$. High-frequency modes (sharp features) die out first; the solution progressively smooths.

Example: The Wave Equation on a String

The same method applies to the wave equation, but with a fundamentally different temporal behavior. Consider a vibrating string fixed at both ends:

Separation Leads to Different Temporal Behavior

Substituting $u(x,t) = X(x) \cdot T(t)$ and separating, the spatial equation is exactly the same as before:

But the temporal equation is second-order:

where $\omega_n = c\sqrt{\lambda_n} = cn\pi/L$ is the natural frequency of mode $n$.

The General Wave Solution

The general solution to the wave equation with zero initial velocity ($g(x) = 0$) is:

with the same Fourier coefficients $b_n = \frac{2}{L}\int_0^L f(x) \sin(n\pi x/L)\,dx$. If the initial velocity is nonzero, we also get sine terms in time with a separate set of coefficients determined by $g(x)$.

Machine Learning Connections

The ideas behind separation of variables are remarkably relevant to modern machine learning. The decomposition of functions into orthogonal modes is the mathematical foundation for several important ML techniques.

1. Fourier Neural Operators (FNO)

The Fourier Neural Operator (Li et al., 2021) directly extends the spectral approach of separation of variables. Instead of solving in physical space, FNO learns solution operators in Fourier (spectral) space:

where $\mathcal{F}$ is the Fourier transform and $R_l$ is a learned weight matrix that acts on the Fourier modes. This is exactly the insight from separation of variables: many PDEs are diagonal in Fourier space, meaning each mode can be processed independently.

2. Spectral Normalization in GANs

Spectral normalization constrains the Lipschitz constant of neural network layers by normalizing their spectral norm (largest singular value). This is an eigenvalue computation, directly analogous to finding eigenvalues in separation of variables. It stabilizes GAN training by controlling how much each "mode" of the discriminator can amplify signals.

3. Principal Component Analysis (PCA)

PCA decomposes data into orthogonal principal components, just as separation of variables decomposes solutions into orthogonal eigenfunctions. The covariance matrix plays the role of the differential operator, eigenvectors play the role of eigenfunctions, and eigenvalues measure the importance of each mode. Keeping only the top modes is like truncating the Fourier series to the most important terms.

4. Diffusion Models

Score-based diffusion models for image generation are built on the same mathematics as the heat equation. The forward diffusion process adds noise progressively — mathematically, this is heat diffusion in pixel space. The Fourier decomposition tells us that high-frequency features (fine details) are destroyed first and low-frequency features (overall structure) last. The reverse (denoising) process reconstructs these features in the opposite order.

Python Implementation

Solving the Heat Equation by Separation of Variables

The following implementation computes the analytical solution using Fourier coefficients and mode superposition — the exact method we derived above.

1import numpy as np
2import matplotlib.pyplot as plt
3
4def separation_of_variables_heat(L=1.0, alpha=0.01, N_modes=20,
5                                  Nx=200, t_values=[0, 0.5, 1, 2, 5]):
6    """
7    Solve the 1D heat equation using separation of variables.
8
9    PDE: du/dt = alpha * d²u/dx²
10    BCs: u(0,t) = u(L,t) = 0
11    IC:  u(x,0) = f(x) (step function)
12
13    Solution: u(x,t) = sum_{n=1}^{N} b_n sin(n*pi*x/L) exp(-alpha*(n*pi/L)^2*t)
14    """
15    x = np.linspace(0, L, Nx)
16
17    # Initial condition: step function
18    def f(x_val):
19        return 1.0 if 0.25 < x_val < 0.75 else 0.0
20
21    # Compute Fourier sine coefficients: b_n = (2/L) * integral f(x) sin(n*pi*x/L) dx
22    def compute_bn(n):
23        integrand = np.array([f(xi) * np.sin(n * np.pi * xi / L) for xi in x])
24        return 2 / L * np.trapezoid(integrand, x)
25
26    # Pre-compute all coefficients
27    coefficients = [compute_bn(n) for n in range(1, N_modes + 1)]
28    print("First 5 Fourier coefficients:")
29    for n, bn in enumerate(coefficients[:5], 1):
30        print(f"  b_{n} = {bn:.6f}")
31
32    # Compute solution at each time
33    fig, axes = plt.subplots(1, len(t_values), figsize=(16, 4))
34
35    for idx, t in enumerate(t_values):
36        u = np.zeros_like(x)
37        for n_idx, bn in enumerate(coefficients):
38            n = n_idx + 1
39            lambda_n = (n * np.pi / L) ** 2
40            # Each mode: b_n * sin(n*pi*x/L) * exp(-alpha*lambda_n*t)
41            u += bn * np.sin(n * np.pi * x / L) * np.exp(-alpha * lambda_n * t)
42
43        axes[idx].fill_between(x, u, alpha=0.3, color='steelblue')
44        axes[idx].plot(x, u, 'b-', linewidth=2)
45        axes[idx].set_title(f't = {t}')
46        axes[idx].set_ylim(-0.3, 1.3)
47        axes[idx].set_xlabel('x')
48        axes[idx].grid(True, alpha=0.3)
49
50    axes[0].set_ylabel('u(x, t)')
51    fig.suptitle(f'Heat Equation Solution ({N_modes} modes)', fontsize=14)
52    plt.tight_layout()
53    plt.show()
54
55    return x, coefficients
56
57# Run the solver
58x, coeffs = separation_of_variables_heat()

Solving the Wave Equation by Separation of Variables

1import numpy as np
2import matplotlib.pyplot as plt
3
4def separation_of_variables_wave(L=1.0, c=1.0, N_modes=30,
5                                  Nx=200, t_values=[0, 0.3, 0.5, 0.8, 1.0]):
6    """
7    Solve the 1D wave equation using separation of variables.
8
9    PDE: d²u/dt² = c² * d²u/dx²
10    BCs: u(0,t) = u(L,t) = 0
11    IC:  u(x,0) = f(x) (triangle), du/dt(x,0) = 0
12
13    Solution: u(x,t) = sum b_n sin(n*pi*x/L) cos(c*n*pi*t/L)
14    """
15    x = np.linspace(0, L, Nx)
16
17    # Initial condition: triangle wave
18    def f(x_val):
19        if x_val < L / 2:
20            return 2 * x_val / L
21        return 2 * (1 - x_val / L)
22
23    # Fourier coefficients
24    def compute_bn(n):
25        integrand = np.array([f(xi) * np.sin(n * np.pi * xi / L) for xi in x])
26        return 2 / L * np.trapezoid(integrand, x)
27
28    coefficients = [compute_bn(n) for n in range(1, N_modes + 1)]
29
30    # Key difference from heat: temporal part is cos, not exp!
31    fig, axes = plt.subplots(1, len(t_values), figsize=(16, 4))
32
33    for idx, t in enumerate(t_values):
34        u = np.zeros_like(x)
35        for n_idx, bn in enumerate(coefficients):
36            n = n_idx + 1
37            omega_n = c * n * np.pi / L  # natural frequency of mode n
38            # Each mode: b_n * sin(n*pi*x/L) * cos(omega_n * t)
39            u += bn * np.sin(n * np.pi * x / L) * np.cos(omega_n * t)
40
41        axes[idx].plot(x, u, 'r-', linewidth=2)
42        axes[idx].fill_between(x, u, alpha=0.2, color='red')
43        axes[idx].set_title(f't = {t}')
44        axes[idx].set_ylim(-1.3, 1.3)
45        axes[idx].set_xlabel('x')
46        axes[idx].grid(True, alpha=0.3)
47
48    axes[0].set_ylabel('u(x, t)')
49    fig.suptitle('Wave Equation: Vibrating String', fontsize=14)
50    plt.tight_layout()
51    plt.show()
52
53    # Show that energy is conserved (unlike heat equation)
54    t_span = np.linspace(0, 3, 300)
55    energy = []
56    for t in t_span:
57        u_vals = np.zeros_like(x)
58        for n_idx, bn in enumerate(coefficients):
59            n = n_idx + 1
60            omega_n = c * n * np.pi / L
61            u_vals += bn * np.sin(n * np.pi * x / L) * np.cos(omega_n * t)
62        # Approximate total energy ~ integral of u^2
63        energy.append(np.trapezoid(u_vals**2, x))
64
65    plt.figure(figsize=(8, 3))
66    plt.plot(t_span, energy, 'g-', linewidth=2)
67    plt.xlabel('Time')
68    plt.ylabel('Energy (L2 norm squared)')
69    plt.title('Wave Equation: Energy is Conserved!')
70    plt.grid(True, alpha=0.3)
71    plt.show()
72
73separation_of_variables_wave()

From Separation of Variables to Spectral Methods

The connection to modern computational methods and ML becomes clear when we implement the spectral approach using FFT:

1import numpy as np
2import matplotlib.pyplot as plt
3
4def spectral_method_pde(N_basis=32, L=2*np.pi, alpha=0.01, dt=0.01, T_final=2.0):
5    """
6    Spectral method for PDEs: the ML-relevant generalization
7    of separation of variables.
8
9    Instead of solving in physical space, we:
10    1. Transform to frequency space (Fourier transform)
11    2. Solve ODEs for each mode independently
12    3. Transform back
13
14    This is the foundation of:
15    - Fourier Neural Operators (FNO)
16    - Spectral normalization in GANs
17    - Frequency-domain learning
18    """
19    x = np.linspace(0, L, N_basis, endpoint=False)
20
21    # Initial condition: sum of two Gaussians
22    u = np.exp(-10 * (x - L/3)**2) + 0.5 * np.exp(-20 * (x - 2*L/3)**2)
23
24    # In spectral space, the heat equation becomes:
25    # d(u_hat_k)/dt = -alpha * k^2 * u_hat_k
26    # This is just an ODE for each mode!
27
28    k = np.fft.fftfreq(N_basis, d=L/N_basis) * 2 * np.pi  # wavenumbers
29
30    # Time evolution in spectral space
31    Nt = int(T_final / dt)
32    history = [u.copy()]
33    times = [0]
34
35    for step in range(Nt):
36        t = (step + 1) * dt
37
38        # Forward FFT: physical -> spectral
39        u_hat = np.fft.fft(u)
40
41        # Exact solution in spectral space (no discretization error!)
42        # Each mode evolves independently: u_hat_k(t+dt) = u_hat_k(t) * exp(-alpha*k^2*dt)
43        u_hat *= np.exp(-alpha * k**2 * dt)
44
45        # Inverse FFT: spectral -> physical
46        u = np.real(np.fft.ifft(u_hat))
47
48        if step % (Nt // 10) == 0:
49            history.append(u.copy())
50            times.append(t)
51
52    # Visualize
53    fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 5))
54
55    for i, (u_snap, t) in enumerate(zip(history, times)):
56        alpha_val = 0.3 + 0.7 * (1 - i / len(history))
57        ax1.plot(x, u_snap, alpha=alpha_val, label=f't={t:.1f}')
58    ax1.set_xlabel('x')
59    ax1.set_ylabel('u(x,t)')
60    ax1.set_title('Spectral Method: Heat Equation Solution')
61    ax1.legend()
62    ax1.grid(True, alpha=0.3)
63
64    # Show spectral decay
65    u_hat_init = np.fft.fft(history[0])
66    u_hat_final = np.fft.fft(history[-1])
67    freqs = np.fft.fftfreq(N_basis, d=L/N_basis)
68    ax2.semilogy(np.abs(freqs[:N_basis//2]),
69                 np.abs(u_hat_init[:N_basis//2]), 'b-o', label='Initial', markersize=3)
70    ax2.semilogy(np.abs(freqs[:N_basis//2]),
71                 np.abs(u_hat_final[:N_basis//2]) + 1e-16, 'r-o', label='Final', markersize=3)
72    ax2.set_xlabel('Frequency |k|')
73    ax2.set_ylabel('|u_hat(k)|')
74    ax2.set_title('Spectral Coefficients: High Frequencies Decay')
75    ax2.legend()
76    ax2.grid(True, alpha=0.3)
77
78    plt.tight_layout()
79    plt.show()
80
81    print("Key insight for ML:")
82    print("  Fourier Neural Operators (FNO) learn solution operators")
83    print("  by parameterizing transformations in spectral space,")
84    print("  just as separation of variables solves PDEs mode-by-mode.")
85
86spectral_method_pde()

Common Pitfalls

Test Your Understanding

Summary

The method of separation of variables is one of the most powerful analytical techniques for solving partial differential equations. By assuming a product form for the solution, we reduce a PDE to a collection of ODEs, each solvable by standard methods.

The Algorithm in Five Steps

1. Assume a product solution: $u(x,t) = X(x) \cdot T(t)$
2. Substitute into the PDE and separate variables to get two sides depending on different variables
3. Set both sides equal to a constant $-\lambda$, obtaining two ODEs
4. Solve the spatial eigenvalue problem using boundary conditions to find eigenvalues $\lambda_n$ and eigenfunctions $X_n(x)$
5. Superpose all modes and determine coefficients from the initial condition using orthogonality

Key Concepts

Key Takeaways

1. Separation of variables reduces a PDE to two or more ODEs, each in a single variable
2. The boundary conditions determine the spatial eigenfunctions, while the initial condition determines the expansion coefficients
3. Orthogonality of eigenfunctions is the key property that makes the method practical — it allows independent computation of each Fourier coefficient
4. The heat equation produces exponential decay of modes (diffusion smooths), while the wave equation produces oscillation (waves vibrate)
5. The same eigenvalue/eigenfunction framework underlies PCA, spectral methods, Fourier Neural Operators, and diffusion models in machine learning
6. The method works for linear PDEs with homogeneous BCs on simple geometries — it does not apply to nonlinear PDEs or irregular domains

Coming Next: In the next chapter, we'll dive deeper into the heat equation, exploring its derivation from physical principles, the heat kernel, and its connections to diffusion in machine learning.