A persistent tension runs through the applied study of partial differential equations (PDEs): the physical world presents us with continuous fields—temperature, pressure, velocity—whose evolution is governed by local relationships, yet the tools we have for understanding those relationships are limited by what we can compute, whether by hand or by machine. The history of the subfield is therefore a story of expanding what counts as a solution, from closed-form formulas to numerical approximations to statistical ensembles, and of the frameworks that made each expansion possible.
The earliest systematic work on PDEs, now gathered under Classical PDE Theory (1740–1900), grew directly out of 18th-century mechanics. The wave equation, the heat equation, and Laplace's equation were the paradigmatic examples. Mathematicians such as d'Alembert, Euler, and Laplace sought explicit formulas—solutions expressed in terms of elementary functions or infinite series—that could be evaluated pointwise. The central pressure was to find a single, deterministic expression that described the entire evolution of a system from its initial and boundary data.
Fourier Analysis and Separation of Variables (1800–1900) emerged as the most powerful technique within this classical tradition. Joseph Fourier's treatment of the heat equation showed that a solution could be represented as a superposition of sine and cosine modes, each evolving independently. This was not merely a computational trick; it introduced the idea that the solution space of a linear PDE could be decomposed into elementary building blocks. Separation of variables became the default method for problems on rectangular, circular, or spherical domains, and it remains a cornerstone of PDE pedagogy today.
Alongside Fourier's approach, the Method of Characteristics (1800–1900) offered a completely different geometric picture. Instead of decomposing the solution into modes, it transformed a first-order PDE into a family of ordinary differential equations along curves called characteristics. For hyperbolic equations such as the wave equation, the method revealed that information propagates along finite-speed trajectories—a physical insight that Fourier analysis alone could not provide. The two frameworks coexisted because they addressed different classes of problems: characteristics for first-order and hyperbolic systems, Fourier series for elliptic and parabolic equations on bounded domains.
Integral Transform Methods (1800–1950) extended the Fourier idea by replacing the discrete series of sines and cosines with continuous transforms—the Fourier transform, the Laplace transform, and later the Hankel and Mellin transforms. These methods were especially effective on unbounded domains, where separation of variables would produce a continuous spectrum rather than a discrete one. They absorbed the logic of Fourier analysis while broadening its reach to problems in diffusion, wave propagation, and potential theory.
A complementary analytical tool arrived with Green's Functions and Integral Equations (1820–1950). George Green's 1828 essay introduced the idea of a point-source response: the solution to a linear PDE with arbitrary forcing could be written as an integral of the forcing against a kernel that depended only on the geometry and the operator. This reframed PDEs as integral equations, shifting attention from the differential operator itself to the properties of the Green's function. The framework coexisted with Fourier methods—Green's functions were often constructed using eigenfunction expansions—but it offered a more direct route to existence proofs and to the treatment of boundary conditions.
By the mid-19th century, the growing catalog of solved problems demanded a unifying classification. The Classification of Second-Order Linear PDEs (1850–1950) sorted equations into elliptic, parabolic, and hyperbolic types based on the sign of the discriminant of their principal part. This was not a mere taxonomy; it predicted the qualitative behavior of solutions—elliptic equations smooth out disturbances instantly, parabolic equations diffuse them, hyperbolic equations propagate them along characteristics. The classification gave applied mathematicians a diagnostic tool: given a new problem, one could first classify it and then choose the appropriate analytical method.
Variational Principles and Weak Solutions (1850–1970) introduced a profound shift in what it meant to solve a PDE. Instead of requiring a twice-differentiable function that satisfies the equation at every point, the variational approach asked for a function that minimizes an energy functional. The Euler–Lagrange equations of the calculus of variations turned out to be the classical PDEs of mechanics and field theory. Later, the concept of a weak solution—a function that satisfies the PDE only when integrated against smooth test functions—allowed mathematicians to handle discontinuous data and irregular geometries. This framework did not replace the classical theory; it enlarged the notion of solution so that existence could be proved under much weaker regularity assumptions. It also laid the conceptual groundwork for the finite element method, which would later discretize the weak form directly.
By the early 20th century, it was clear that most PDEs of practical interest—especially nonlinear ones—could not be solved by any of the classical analytical methods. Asymptotic and Perturbation Methods (1900–Present) addressed this gap by constructing approximate solutions in limiting regimes. When a small parameter appeared in the equation, one could expand the solution as a power series in that parameter, matching terms order by order. Boundary-layer theory in fluid dynamics, the WKB approximation in quantum mechanics, and multiple-scale analysis in nonlinear oscillations all grew from this framework. Asymptotic methods did not compete with numerical methods; they complemented them by providing analytical insight into the structure of solutions—the location of shocks, the thickness of layers, the frequency of oscillations—that pure computation could not easily reveal.
Finite Difference Methods (1920–Present) were the first systematic numerical framework for PDEs. The idea was straightforward: replace derivatives with difference quotients on a grid, turning the PDE into a system of algebraic equations. Early work by Courant, Friedrichs, and Lewy in 1928 established the CFL condition, which linked the grid spacing to the time step for stability. Finite difference methods were simple to implement and worked well on rectangular domains, but they struggled with complex geometries and with conservation laws where naive differencing could produce spurious oscillations. They coexisted with asymptotic methods throughout the mid-20th century, each handling the problems the other could not.
Finite Element Methods (1950–Present) emerged from structural engineering and were placed on a rigorous mathematical foundation by the 1960s. Instead of approximating derivatives on a grid, the finite element method partitions the domain into small elements (triangles, quadrilaterals, tetrahedra) and represents the solution as a piecewise polynomial that satisfies the weak form of the equation. This framework absorbed the variational principle directly: the weak form became the starting point, and the finite-dimensional subspace of test functions became the computational approximation. Finite elements excelled where finite differences struggled—complex geometries, heterogeneous materials, and problems with natural boundary conditions. Over time, the two frameworks have coexisted with a rough division of labor: finite differences remain popular for simple geometries and time-dependent problems, while finite elements dominate structural mechanics, electromagnetics, and multiphysics simulations.
Spectral and Pseudospectral Methods (1970–Present) took the Fourier idea to its computational extreme. Instead of using local basis functions (as in finite elements) or local differences (as in finite differences), spectral methods represent the solution as a global expansion in smooth basis functions—Fourier series on periodic domains, Chebyshev polynomials on bounded intervals. When the solution is smooth, spectral methods achieve exponential convergence with far fewer degrees of freedom than any local method. They are the method of choice for problems in turbulence, weather prediction, and wave propagation where accuracy is paramount. Their limitation is that they degrade sharply in the presence of discontinuities or complex geometry, so they coexist with finite elements and finite differences rather than replacing them.
Soliton Theory and Integrable Systems (1960–Present) revealed that some nonlinear PDEs possess an extraordinary hidden structure: they can be solved exactly by a method called inverse scattering, which generalizes the Fourier transform to nonlinear problems. The Korteweg–de Vries equation, the nonlinear Schrödinger equation, and the sine-Gordon equation all turned out to be integrable, supporting solitary waves (solitons) that pass through each other without change of shape. This framework revived the classical ideal of exact solutions, but for a very special class of nonlinear equations. It coexists with numerical and asymptotic methods because most nonlinear PDEs are not integrable; soliton theory provides deep insight into the exceptional cases that are.
Stochastic PDEs (1970–Present) introduced randomness into the equation itself. Instead of a deterministic forcing or initial condition, the PDE includes a noise term—typically white noise in space and time—so that the solution becomes a random field. This framework was driven by problems in statistical physics (fluctuating interfaces, reaction-diffusion systems), finance (interest rate models), and fluid dynamics (turbulent transport). Stochastic PDEs required new analytical tools—the Itô calculus extended to infinite dimensions, the theory of martingale solutions—and new numerical methods that could handle the roughness of the noise. They coexist with deterministic PDE theory as a parallel track: the deterministic frameworks handle the mean behavior, while stochastic PDEs capture the fluctuations around it.
Data-Driven PDE Methods (2010–Present) represent the most recent expansion of what counts as a solution. When the governing equation is unknown or too complex to write down, one can learn an approximate PDE from observational data using neural networks, Gaussian processes, or sparse regression. The physics-informed neural network (PINN) framework, for example, incorporates the residual of a known PDE as a penalty term in the loss function, allowing the network to approximate solutions without a grid. Other approaches discover the PDE itself from data, identifying the terms and coefficients that best explain the observed evolution. Data-driven methods do not replace the classical or numerical frameworks; they extend them to regimes where the equation is uncertain, the geometry is irregular, or the data are sparse. They are currently the most active area of methodological development in the subfield.
The leading frameworks today—Asymptotic and Perturbation Methods, Finite Difference Methods, Finite Element Methods, Spectral and Pseudospectral Methods, Soliton Theory, Stochastic PDEs, and Data-Driven PDE Methods—coexist in a productive division of labor. They agree on the central importance of the classification of PDEs and on the variational principle as a unifying concept. They disagree on what constitutes a satisfactory solution: asymptotic analysts prize structural insight, numerical analysts prioritize accuracy and efficiency, stochastic analysts emphasize ensemble behavior, and data-driven practitioners value flexibility and automation. The classical analytical frameworks (Fourier analysis, Green's functions, separation of variables) remain essential for teaching, for benchmarking numerical methods, and for problems where exact solutions are still attainable. The tension between continuous ideals and computational realities that opened the subfield's history has not been resolved; it has instead become the engine that drives the development of new frameworks, each expanding the range of PDEs that applied mathematicians can understand and simulate.