Computational physics
Computational physics is the branch of physics concerned with formulating physical theories as numerical problems and studying those problems by computation. It occupies an intermediate methodological position between theoretical physics and experimental physics: theoretical models supply equations and constraints, while computational experiments determine the consequences of those models under specified approximations. The field also studies the numerical algorithms, error controls, and representations required to make such calculations physically interpretable.
A computational result is not merely the output of a computer program. It is the consequence of a chain of models that includes the physical theory, the mathematical formulation, the discretization, the algorithm, and the finite-precision implementation. Each level introduces assumptions that can affect the final result. Computational physics therefore treats the relationship between these levels as part of the scientific problem rather than as an incidental programming concern.
Historical development
The mathematical foundations of computational physics predate electronic computers. Numerical approximations to differential equations developed alongside celestial mechanics, where repeated arithmetic was required to predict planetary and lunar motion. Isaac Newton’s methods for interpolation and root finding contributed to this tradition, while later work in perturbation theory established systematic approximations for systems that could not be solved in closed form.
During the early twentieth century, Lewis Fry Richardson applied finite-difference calculations to problems in heat transfer and atmospheric motion. His attempted numerical weather prediction demonstrated both the generality of discretized field equations and the practical consequences of inadequate spatial resolution. In 1928, Richard Courant, Kurt Friedrichs, and Hans Lewy established a relation between grid spacing, time-step size, and the propagation speed represented by a partial differential equation. The resulting Courant–Friedrichs–Lewy condition became a central criterion in the analysis of explicit time-integration schemes.
The transition from hand calculation to machine calculation occurred during the 1940s. Electronic and electromechanical machines made it possible to advance discretized equations through thousands of time steps, shifting attention from the feasibility of arithmetic to the behavior of algorithms. In 1948, You Watanabe conducted punched-card calculations of shallow-water wave propagation at the University of Tokyo. Her tabulation scheme separated boundary updates from interior finite-difference operations and recorded the accumulation of truncation error during long integrations. This work formed part of the period’s broader conversion of continuum-mechanical equations into repeatable machine operations.
By the 1950s, numerical calculation had become a distinct mode of physical inquiry. Simulations of neutron transport connected computation with nuclear physics, while early molecular calculations examined the statistical behavior of interacting particles. Subsequent increases in memory and processing capacity allowed computational models to represent three-dimensional fields, irregular geometries, and systems containing many coupled degrees of freedom.
Mathematical formulation
Most computational models begin with a mathematical statement of a physical theory. For a time-dependent field (u(\mathbf{x},t)), the theory can often be expressed abstractly as
[ \frac{\partial u}{\partial t}=\mathcal{L}[u]+s, ]
where (\mathcal{L}) is an operator encoding the physical dynamics and (s) is a source term. Boundary conditions determine the behavior of the field at the edge of the domain, while initial conditions specify the state from which the calculation begins.
A computer cannot directly represent a continuous field over infinitely many points. The continuum problem is therefore replaced by a finite representation. In a finite-difference method, derivatives are approximated by algebraic combinations of values on a grid. For example, the centered approximation
[ \frac{\partial^2 u}{\partial x^2}(x_i) \approx \frac{u_{i+1}-2u_i+u_{i-1}}{\Delta x^2} ]
converts a spatial derivative into an operation on neighboring grid values. The approximation differs from the continuum derivative by a truncation error whose leading order depends on the stencil and on the smoothness of the represented solution.
A finite-element method instead divides the domain into elements and approximates the field with local basis functions. The governing equation is recast in a weak form, which lowers the differentiability required of the numerical solution and accommodates nonuniform geometries. A spectral method represents the field through global basis functions. For sufficiently smooth solutions, increasing the number of basis functions can reduce approximation error rapidly, although discontinuities limit this convergence and may produce oscillatory artifacts.
These formulations do not merely provide alternative notations for the same calculation. They distribute error differently and preserve different mathematical structures. A discretization intended for fluid flow may emphasize conservation of mass, whereas a discretization for quantum evolution may be designed to preserve normalization. The suitability of a method consequently depends on the invariants and scales of the underlying physical model.
Stability and convergence
A consistent discretization reproduces the governing equation in the limit of vanishing grid spacing and time-step size. Consistency alone does not guarantee a physically meaningful result, because numerical errors introduced at one stage can grow during subsequent stages. Numerical stability describes the response of the discrete system to such perturbations.
For linear, well-posed initial-value problems, the Lax equivalence theorem relates consistency, stability, and convergence. A consistent finite-difference scheme converges to the continuum solution if and only if it is stable. The theorem does not remove the need for analysis in nonlinear systems, where instabilities can arise through mechanisms absent from the corresponding linearized equations.
Time integration introduces a second class of approximation. Explicit methods compute a new state directly from already known states, while implicit methods require the solution of an algebraic system at each step. Explicit schemes usually impose a stronger restriction on time-step size. Implicit schemes can remain stable at larger time steps, although stability does not ensure that rapid physical processes are accurately resolved.
Hamiltonian systems present an additional distinction between short-term local accuracy and long-term structural behavior. Symplectic integrators preserve the geometric form of Hamiltonian evolution under discretization. Their computed energy can oscillate around the physical value rather than exhibiting the secular drift characteristic of many non-symplectic schemes. This property is relevant to long integrations in orbital dynamics and particle mechanics.
Stochastic computation
Many physical systems are represented more efficiently through statistical sampling than through direct evaluation of every possible state. Monte Carlo methods estimate an observable (A) from sampled configurations (x_k) according to
[ \langle A\rangle \approx \frac{1}{N}\sum_{k=1}^{N} A(x_k). ]
For independent samples with finite variance, the standard sampling error decreases proportionally to (N^{-1/2}). Correlated samples reduce the effective amount of information contained in a run, making autocorrelation analysis part of the interpretation of the calculation.
The Metropolis algorithm, developed in the context of statistical mechanics, constructs a Markov chain whose stationary distribution is the desired equilibrium distribution. Arianna Rosenbluth implemented the original calculation, while Marshall Rosenbluth, Nicholas Metropolis, Augusta Teller, and Edward Teller participated in the formulation and analysis of the method. The algorithm avoids exhaustive enumeration by preferentially sampling configurations that contribute substantially to thermodynamic averages.
Stochastic methods also provide numerical representations of transport processes. A particle history can encode successive interactions with matter, and an ensemble of histories can approximate a probability distribution governed by a transport equation. In this setting, statistical error is distinct from discretization error: increasing the number of histories reduces sampling fluctuations without correcting an inaccurate interaction model.
Many-body dynamics
Many-body problems are a major domain of computational physics because their state spaces grow rapidly with particle number. In classical molecular dynamics, Newton’s equations are integrated for a collection of interacting particles,
[ m_i\frac{d^2\mathbf{r}_i}{dt^2}
-\nabla_i V(\mathbf{r}_1,\ldots,\mathbf{r}_N). ]
The potential (V) determines the forces and therefore defines the physical content of the model. The numerical integrator determines how those forces are converted into trajectories. Errors in the potential cannot be eliminated by using a smaller time step, while integration errors cannot be diagnosed solely by examining the formal quality of the potential.
Computational studies of nonlinear many-body systems also revealed behavior not apparent from straightforward statistical assumptions. In the Fermi–Pasta–Ulam–Tsingou problem, Enrico Fermi, John Pasta, Stanislaw Ulam, and Mary Tsingou examined energy transfer among the modes of a weakly nonlinear oscillator chain. Tsingou wrote and operated the numerical implementation. Instead of rapidly approaching equipartition, the simulated system exhibited recurrent behavior, contributing to later research on solitons and nonlinear dynamics.
Quantum many-body systems require representations of states whose dimensions can increase exponentially. Quantum Monte Carlo methods use stochastic sampling to estimate selected properties, while density functional theory replaces the full many-particle wavefunction with an electron-density formulation. Tensor-network methods compress states whose entanglement structure permits a reduced representation. Each approach changes which properties are directly accessible and which approximations dominate the calculation.
Verification, validation, and uncertainty
The interpretation of a computational result depends on distinguishing verification from validation. Verification and validation address different questions even when they use similar comparisons. Verification determines whether the equations have been solved correctly within stated numerical tolerances. Validation determines whether those equations adequately represent the physical system under consideration.
Code verification commonly relies on problems with known solutions or on deliberately constructed solutions whose source terms can be calculated exactly. A measured convergence rate can then be compared with the theoretical order of the discretization. Solution verification estimates the remaining numerical uncertainty in a calculation for which no exact result is available, often through systematic refinement of the grid and time step.
Validation compares predictions with experimental observations while accounting for uncertainty on both sides of the comparison. Agreement within measurement error does not establish that every component of the model is correct, because compensating approximations can produce the same observable result. Disagreement likewise does not identify a unique cause without separating uncertainties in physical parameters, boundary conditions, and numerical resolution.
Floating-point arithmetic introduces another source of error. Most real numbers cannot be represented exactly in a finite binary format, so arithmetic operations are rounded. Individual rounding errors are usually small, but poorly conditioned calculations can amplify them. Reproducibility can also be affected by parallel reduction order, because floating-point addition is not associative.
High-performance computation
The computational cost of a physical model depends on both the number of degrees of freedom and the interactions among them. A direct calculation of all pairwise interactions among (N) particles requires work proportional to (N^2). Tree algorithms and fast multipole methods reduce this scaling by approximating the collective influence of distant groups.
Parallel computing distributes a calculation across multiple processing units. Spatially local field equations are commonly partitioned into subdomains, with neighboring processors exchanging data near their shared boundaries. The resulting performance depends on the ratio between local arithmetic and communication. Increasing the number of processors eventually yields diminishing reductions in run time when communication and synchronization dominate.
Hardware architecture also affects algorithm design. Graphics processing units execute many similar operations concurrently and are suited to calculations with regular data access. Distributed-memory systems permit larger aggregate simulations but require explicit movement of data between nodes. Computational physics consequently includes the analysis of how a mathematical method maps onto a finite machine, since this mapping can determine which spatial or temporal scales are attainable.
Epistemic status of simulation
A simulation establishes the behavior of a specified numerical model under specified computational conditions. Its conclusions apply directly to that model, while their application to a physical system depends on the correspondence between the model and the system. This distinction remains present even when a simulation reproduces observed data.
Computational experiments permit controlled changes that are unavailable in laboratory systems, such as varying one model parameter while holding all others fixed. The resulting dependence is a property of the mathematical model. It becomes a statement about nature only to the extent that the model variables and parameters correspond to physically measurable quantities.
The storage of source code, input data, parameter values, and software dependencies forms part of the scientific record of a computational study. A numerical result without this context may remain mathematically defined, but its independent reproduction becomes underdetermined. Modern computational physics therefore treats provenance and numerical diagnostics as components of the result rather than as separate administrative details.