# Computational physics **Computational physics** is the practice of solving physical problems by stepping them forward on a machine rather than solving them in closed form. Almost no equation of interest in [[Physics|physics]] has an exact solution — the pendulum with a real amplitude, three gravitating bodies, a lattice of interacting spins, a turbulent flow — and the discipline exists because a discretised, approximate answer computed in finite time is worth more than an exact answer nobody can write down. Its central and most uncomfortable fact is that the approximation has two independent errors moving in opposite directions: cut the timestep and truncation error falls, but the number of arithmetic operations rises and their rounding accumulates. In the microsim below the reader has one slider, the number of timesteps N_t on a logarithmic axis from 50 to 100,000, and one toggle marked *float*. Four integrators — Euler, Euler–Cromer, [[Verlet_integration|Verlet]] and fourth-order [[Runge–Kutta_methods|Runge–Kutta]] — run the same [[Pendulum|pendulum]], `theta'' = −(g/l)·sin(theta)`, with g/l = 10, theta0 = 0.2 rad and v0 = 0 over t ∈ [0, 6], against the small-angle exact curve `x = x0·cos(w·t) + (v0/w)·sin(w·t)`.[^anag-motion] Three panels answer: theta(t) laid over the exact curve, log10 of the pointwise error for each method, and the final error plotted against N_t, whose slopes are the whole lesson — one decade of error per decade of N_t for the two Euler methods, two for Verlet, four for RK4 — until every curve flattens onto a roundoff floor near N_t ≈ 50,000.[^anag-motion][^manual04] The *float* toggle pushes the state through single precision and the floor jumps forward to N_t ≈ 18,000, reproducing the failure the book prints.[^anag-motion][^spec-p46] On the [[Physics]] flagship this article serves *Current research* in Part III — Research: it is the page that says how physics is actually computed today, and the open questions listed below are the ones being attacked with these tools. ## Overview Computational physics sits between analytic theory and experiment and borrows the vocabulary of both. The analyst writes a differential equation; the computational physicist replaces the derivative by a difference, choosing a step and a stencil, and then has to prove that the sequence of numbers produced has anything to do with the solution. For the pendulum the first-order form is `dx/dt = v`, `dv/dt = a(t, x, v)`, and the simplest possible rule — advance both by the current rate over a step `dt = (t_f − t_i)/(N_t − 1)` — is the Euler method, whose local error is O(dt²) and whose global error is O(dt).[^anag-motion] One line of the sim's code differs between Euler and Euler–Cromer: the latter uses the *updated* velocity to advance the position. The two have the same order, and on the error-against-N_t plot they lie almost on top of each other, yet over a long conservative run their behaviour is not the same at all — which is the first hint that order alone is the wrong figure of merit.[^manual02] The numbers are checkable, and in this field they must be. Anagnostopoulos prints the first rows of his `euler.dat` for x0 = 0.2, v0 = 0, t_f = 6, N_t = 1000: `0 0.20000 0`, then `0.00600 0.20000 −0.01193`, `0.01201 0.19992 −0.02386`, `0.01801 0.19978 −0.03579`.[^anag-motion] Those rows reproduce exactly when the step is 6/999 = 0.006006 and the acceleration carries the full sin(theta) rather than its small-angle approximation, and the listing truncates rather than rounds, so a comparison must allow 10⁻⁵.[^anag-motion][^derived-cp] With omega² = 10 the period is T = 2π/√10 = 1.987 s, so the six-second window holds a little over three swings — long enough for an unstable method to show itself and short enough to draw.[^anag-motion][^derived-cp] ### Status in physics Computational work is now a third mode of practice rather than a service function, with its own failure modes, its own literature and its own idea of a reproducible result. It is also where the discipline's open problems are attacked. No settled microscopic theory explains [[High-temperature_superconductivity|high-temperature superconductivity]] in the cuprates; the identity of [[Dark_matter|dark matter]] remains unresolved after decades of direct and indirect search; [[Quantum_gravity|quantum gravity]] has no experimentally tested formulation; and [[Turbulence|turbulence]] — a classical problem, in a classical equation, written down in 1845 — still has no closure that predicts a flow from first principles.[^openprob] Each of the four is a computational problem in practice: lattice simulations of correlated electrons, N-body and structure-formation runs, discrete approaches to spacetime, and direct numerical simulation of the Navier–Stokes equations. All four appear on the physics community's standing list of unsolved problems, of which [[List_of_unsolved_problems_in_physics|the published catalogue]] is the convenient index. ## Challenges in computational physics The characteristic difficulty is that making a calculation more accurate eventually makes it worse. Every floating-point operation rounds, and the size of that rounding is set by machine epsilon: 2.2×10⁻¹⁶ in double precision and 1.2×10⁻⁷ in single.[^anag-motion] Truncation error falls as a power of dt, so it falls as N_t rises; accumulated roundoff rises with the number of operations, so it rises with N_t. Where the two cross, the total error stops improving and the curve goes flat. In the sim that floor sits near N_t ≈ 50,000 in double precision for Verlet, and the *float* toggle moves it to roughly 18,000 — the reader can watch a fourth-order method lose to a second-order one simply because it has run out of digits.[^anag-motion][^manual04] Stability is a separate and sharper failure. On the linear test equation `y' = lambda·y`, each explicit method multiplies the solution by an amplification factor per step, `Q_FE = 1 + z`, `Q_TR = (1 + z/2)/(1 − z/2)`, `Q_RK4 = 1 + z + z²/2 + z³/6 + z⁴/24` with z = lambda·dt, and the method is stable only while |Q| ≤ 1.[^vuik] That single inequality sets hard step bounds: forward Euler and modified Euler need dt ≤ −2/lambda, RK4 needs dt ≤ −2.8/lambda, and the implicit backward-Euler and trapezoidal rules have no bound at all.[^vuik] For the stiff problem `y' = −100(y − cos t) − sin t`, forward Euler is limited to dt = 0.02 for the whole run by a transient that is physically dead after a twentieth of a second — the definition of stiffness.[^vuik] On the imaginary axis, where an undamped oscillator lives, forward Euler is unstable at *every* step size, and RK4 survives only while |lambda·dt| ≤ 2√2 = 2.828.[^vuik][^derived-cp] Nonlinearity adds a third. The stability bound is local, derived from a local value of lambda, so a method that is stable on a small swing can fail on a large one: at theta0 = 3.0, close to inverted, Euler is unstable for N_t below about 1000 on the same six-second window.[^anag-motion] Beyond that lies genuine [[Chaos_theory|chaos]], where two runs differing in the last bit separate exponentially and only statistical statements survive. ## Methods and algorithms Four stepping rules cover most of what the sim shows. Euler and Euler–Cromer are first order. Euler–Verlet, `theta_{n+1} = 2·theta_n − theta_{n−1} + alpha_n·dt²`, is second order, needs a seed step `theta_0 = theta_1 − omega_1·dt + (1/2)·alpha_1·dt²`, and requires a one-sided difference for the final velocity.[^anag-motion] Fourth-order Runge–Kutta evaluates the acceleration four times per step and is globally O(dt⁴).[^anag-motion] Those orders are exactly the slopes on the sim's third panel: one decade of error removed per decade of N_t for both Euler variants, two for Verlet and four for RK4.[^anag-motion] The cost side is where order pays. To reach an error of 10⁻⁴ on the unit interval, forward Euler needs dt = 10⁻⁴ and about 10,000 function evaluations, modified Euler needs dt = 10⁻² and 200, and RK4 needs dt = 10⁻¹ and 40 — a factor of 250 between the first and the last for the same answer.[^vuik] This is why nobody writes a production solver in Euler, and why the sim is honest in showing Euler at all: it is the baseline against which the others are measured. Order is nevertheless not the whole story, and the sim's fourth lesson is structural. For a conservative system, a [[Symplectic_integrator|symplectic]] method such as Verlet does not conserve the true energy either, but it conserves a nearby shadow energy, so E oscillates with an O(dt²) amplitude instead of drifting. RK4 has no such property: on a long orbit it drifts slowly and monotonically, and an adaptive stepper drifts too — tightening a relative tolerance from 10⁻³ to 10⁻⁵ slows the drift without removing it.[^downey-orbit][^manual02] Forward Euler on the harmonic oscillator is worse than drifting: its amplification factor multiplies the energy by 1 + (omega·dt)² every step, so the orbit spirals outward by construction.[^vuik][^derived-cp] For an undamped run, choosing Verlet over RK4 matters more than choosing fourth order over second. ## Divisions The field divides by the equation being discretised rather than by the machine. Molecular and [[Materials_science|materials]] work runs [[Molecular_dynamics|molecular dynamics]], where the same Verlet step that drives the sim is applied to 10⁶ atoms and where energy conservation over a nanosecond is the acceptance test. [[Computational_chemistry|Computational chemistry]] and [[Density_functional_theory|density functional theory]] solve the electronic structure problem instead, turning a many-body [[Schrödinger_equation|Schrödinger equation]] into a self-consistent field. [[Statistical_mechanics|Statistical mechanics]] is computed stochastically: the [[Ising_model|two-dimensional Ising model]] under [[Metropolis–Hastings_algorithm|Metropolis]] dynamics is the standard laboratory, and a well-written implementation caches the small set of possible acceptance ratios in a lookup table rather than calling the exponential per spin.[^anag-ising] [[Fluid_dynamics|Fluid]] and plasma work discretises the continuum equations in space as well as time, as does [[Magnetohydrodynamics|magnetohydrodynamics]]; [[Astrophysics|astrophysical]] N-body and [[Condensed_matter_physics|condensed-matter]] lattice simulations sit at the two extremes of scale. Across all of them the [[Monte_Carlo_method|Monte Carlo method]] appears wherever a high-dimensional integral has to be sampled rather than evaluated.[^anag-ising] ## Applications The applications that matter are the ones where no experiment is available. An [[Orbit|orbit]] integrated over a billion years, a [[Nuclear_fusion|fusion]] plasma at conditions no laboratory has yet reached, a protein folding on a microsecond that instruments cannot resolve — each is a computation standing in for a measurement, and each inherits the error budget above. Downey's textbook example is deliberately small and deliberately damning: a two-body gravitational orbit with `KE = m·v²/2` and `PE = −G·m1·m2/r`, integrated by a general-purpose adaptive solver, does not conserve energy, and the drift is visible long before anything else in the run looks wrong.[^downey-orbit] The lesson generalises. Any result from a long conservative integration must be reported with the energy history beside it, because a plausible trajectory computed by an energy-leaking method is the most persuasive wrong answer in the field. Where the underlying system is a [[Cellular_automaton|cellular automaton]] or an [[Agent-based_model|agent-based model]] rather than a differential equation, the analogous check is a conserved count. ## Software Practically, computational physics is written in a compiled language for the inner loop and a scripted one for everything around it, and the textbooks reflect the split: Anagnostopoulos develops the integrators in C++ with full source listings and plotted output, while Downey teaches the same material through a high-level environment's built-in solvers.[^anag-motion][^downey-ode] Both routes carry traps that are software, not physics. The RK4 step in the C++ text should be coded from the book's `RKSTEP` listing rather than from the typeset equation, which was damaged in extraction.[^anag-motion][^manual04] A general-purpose adaptive solver hides its step-size logic, which is convenient until the energy drift has to be explained.[^downey-orbit] In a browser sim the corresponding trap is the storage type: the integrators operate on typed arrays, and a `Float32Array` state silently reproduces the book's single-precision failure near N_t ≈ 18,000 — which is why the Wikitube build fixes `Float64Array` and exposes single precision only as a deliberate toggle.[^manual04] ## See also - [[Verlet_integration]] — the second-order symplectic step used by most of the cluster's sims - [[Runge–Kutta_methods]] - [[Symplectic_integrator]] - [[Molecular_dynamics]] - [[Monte_Carlo_method]] - [[Ising_model]] - [[Ordinary_differential_equation]] - [[List_of_unsolved_problems_in_physics]] ## References [^anag-motion]: Anagnostopoulos, Konstantinos (2016). *Computational Physics: A Practical Introduction to Computational Physics and Scientific Computing (using C++)*, 2nd ed. Chapter 6, *Motion of a Particle*: the first-order form `dx/dt = v`, `dv/dt = a(t, x, v)`, p. 202–203; the pendulum `theta'' = −(g/l)·sin(theta)`, p. 203; Euler's local O(Δt²) and global O(Δt) error with `dt = (t_f − t_i)/(N_t − 1)`, p. 204; Euler–Cromer using the updated velocity, p. 205; the Euler–Verlet recurrence, its half-step seed and its one-sided final velocity, pp. 204–206; RK2 global O(h²) and RK4 global O(h⁴), pp. 217–219; the `RKSTEP` listing to be coded in place of the typeset Eq. 4.22, pp. 220, 223; the harmonic reference solution and energy, pp. 224–225; machine epsilon 2.2×10⁻¹⁶ (double) and 1.2×10⁻⁷ (float) and the single-precision Verlet floor near N_t ≈ 18,000, pp. 206–209; the `euler.dat` opening rows for x0 = 0.2, v0 = 0, t_f = 6, N_t = 1000, p. 213; the error falling 10×, 10×, 100× and 10⁴× per decade of N_t until roundoff near N_t = 50,000, pp. 224–225; ω² = 10 giving T ≈ 1.987, p. 205; Euler unstable at θ0 = 3.0 for N_t ≲ 1000, p. 208. [^anag-ising]: Anagnostopoulos, Konstantinos (2016). *Computational Physics*, 2nd ed. Chapter 14, *Monte Carlo Simulations*, pp. 500–519, and Chapter 15, *Simulation of the d = 2 Ising Model*, pp. 520–593: the Metropolis implementation and its acceptance lookup table. [^vuik]: Vuik, Kees; Vermolen, Fred; van Gijzen, Martin (2023). *Numerical Methods for Ordinary Differential Equations*. Chapter 6, *Numerical time integration of ODEs*: the test equation `y' = lambda·y`, pp. 78–79; forward Euler, backward Euler, trapezoidal and modified Euler, pp. 76–77, and RK4, pp. 85–86; the amplification factors and the |Q| ≤ 1 stability criterion, pp. 79–80, 86; the step bounds Δt ≤ −2/λ for FE and ME and Δt ≤ −2.8/λ for RK4, with BE and TR unconditionally stable, pp. 80–81, 86–87; method orders 1, 1, 2, 2 and 4, pp. 83, 87; behaviour on the imaginary axis, where FE and ME are always unstable and RK4 needs |λΔt| ≤ 2.8, pp. 89, 99–100; the local-λ bound for nonlinear problems, p. 82; the stiff example `y' = −100(y − cos t) − sin t` with |Q_BE| = 1/21 and |Q_TR| = 9/11 at Δt = 0.2 against an FE limit of Δt = 0.02, pp. 102–103; and the cost to reach 10⁻⁴ on [0, 1] — FE 10,000 evaluations, ME 200, RK4 40 — pp. 85, 91. [^downey-orbit]: Downey, Allen (2021). *Physical Modeling in MATLAB*, version 4.0. Chapter 15, *Springs and Things*, pp. 155–160: the inverse-square force with the Sun held fixed and the energies `KE = m·v²/2`, `PE = −G·m1·m2/r`, p. 159; and the finding that an adaptive `ode45` orbit does not conserve energy, that tightening the relative tolerance from 10⁻³ to 10⁻⁵ slows but does not stop the drift, and that `ode23` "works surprisingly well", p. 160. [^downey-ode]: Downey, Allen (2021). *Physical Modeling in MATLAB*, version 4.0. Chapter 10, *Ordinary Differential Equations*, pp. 99–112: Euler's method presented as the assumption that the rate is constant across a step, p. 101, with the yeast-growth test in which Euler at Δt = 0.1 and the adaptive solver differ by under 1 % at t = 4, pp. 102–104. [^openprob]: The four open problems named here are cited to their standard reference works rather than to a Portal Book, none of which covers them: Anderson, P. W. (1997). *The Theory of Superconductivity in the High-Tc Cuprates*, Princeton University Press (page to pin), for the absence of a settled microscopic theory of cuprate superconductivity; Particle Data Group, *Review of Particle Physics*, review article on dark matter (page to pin); Rovelli, Carlo (2004). *Quantum Gravity*, Cambridge University Press (page to pin), for the state of the problem; and Frisch, Uriel (1995). *Turbulence: The Legacy of A. N. Kolmogorov*, Cambridge University Press (page to pin), for the closure problem in the Navier–Stokes equations. [^manual04]: Wikitube MicroSim Guide, sub-manual 04 *Atomic, Quantum, Statistical and Electromagnetic Physics*, §6.1 "Integrator error: Euler, Euler–Cromer, Verlet, RK4": the microsim concept and its single control (N_t, log, 50–100,000, default 1000) with a `float` toggle; the three panels; the HUD line `dt = (tf - ti)/(Nt - 1); err ~ dt^p`; and the build pitfalls, including the requirement to use `Float64Array` because a `Float32Array` state quietly reproduces the book's single-precision failure, and the instruction to seed Verlet with the half-step or lose its second-order accuracy. [^manual02]: Wikitube MicroSim Guide, sub-manual 02 *Mechanics*, §5.1 "Stability regions and stiffness" and §5.2 "Energy drift: orbits and symplectic stepping": the stiff race and its rasterised stability regions; the statement that a second-order symplectic step beats fourth-order RK4 on energy over long runs, so that the choice of `verlet` over `rk4` matters more than the order; and the note that Verlet conserves a shadow energy, with E oscillating at O(Δt²) amplitude rather than staying exactly constant. [^spec-p46]: Matter & Energy Cluster contract, `_registry/plans/PHYSICS_SECTIONS.md` row P46: the sim concept for this page — four methods run together on `theta'' = −(g/l)·sin(theta)` with ω² = 10, theta0 = 0.2 and t ∈ [0, 6]; N_t as the single control on a log axis from 50 to 100,000; x(t) against the exact curve, log10 of the per-method error, and the error-against-N_t line with slopes 1, 1, 2 and 4 down to the roundoff floor near N_t ≈ 50,000; the `float` toggle reproducing the early floor at N_t ≈ 18,000; and the instruction that the open problems be carried as prose with links rather than as sim state. [^derived-cp]: Computed for this article from the equations and constants cited above: the step `dt = 6/999 = 0.0060060` s that reproduces the printed `euler.dat` times 0.00600, 0.01201, 0.01801; the period `T = 2π/√10 = 1.98692` s at ω² = 10, giving 3.02 periods in the six-second window; the RK4 imaginary-axis stability bound |λΔt| ≤ 2√2 = 2.8284, which the source rounds to 2.8; the forward-Euler energy factor on `y'' = −y`, `|Q|² = 1 + (ωΔt)²` per step, from `Q_FE = 1 + z` with z = ±iωΔt; and the evaluation-count ratio 10,000 ÷ 40 = 250 between forward Euler and RK4 at a target error of 10⁻⁴. ## Further reading - Anagnostopoulos, Konstantinos (2016). *Computational Physics: A Practical Introduction to Computational Physics and Scientific Computing (using C++)*, 2nd ed. CC BY-SA. Chapters 6 and 15 are the source of the integrator study and the Ising implementation used here. - Vuik, Kees; Vermolen, Fred; van Gijzen, Martin (2023). *Numerical Methods for Ordinary Differential Equations*. CC BY. Chapter 6 is the stability-region treatment. - Downey, Allen (2021). *Physical Modeling in MATLAB*, version 4.0. CC BY-NC. Chapters 10 and 15 for the adaptive-solver view and the orbital energy-drift example. ## External links - [*Computational Physics* (Anagnostopoulos, 2016)](https://open.umn.edu/opentextbooks/textbooks/computational-physics-a-practical-introduction-to-computational-physics-and-scientific-computing-using-c) — the second edition, with the C++ listings behind this page - [*Physical Modeling in MATLAB* (Downey, 2021)](https://open.umn.edu/opentextbooks/textbooks/physical-modeling-in-matlab) — the orbit and ODE chapters - The Wikipedia pair's *External links* section lists the discipline's societies, journals and code repositories. <!-- MATTERSIM:BEGIN g24 — Matter & Energy Cluster microsim (framework build, specs/sims/Computational_physics.json); do not hand-edit inside --> **Microsim — three.js (Wikitube framework), pending deploy:** *Computational physics* will play here once `https://wikitube-3d-microsims.netlify.app/matter/Computational_physics.html` is live. <!-- pending: <div class="wt-sim" data-src="https://wikitube-3d-microsims.netlify.app/matter/Computational_physics.html" data-title="Computational physics"></div> --> <!-- MATTERSIM:END --> ## Wikipedia : Wikitube **Strict pair:** [Wikipedia](https://en.wikipedia.org/wiki/Computational_physics) : [Wikitube](https://en.wikitube.io/wiki/Computational_physics) · pinned revision [1323914937](https://en.wikipedia.org/w/index.php?oldid=1323914937) · 2026-09-11 ## Previous hub tags Hubs: `Life_Physics`. Portals: [[PORTAL_Physics]]. --- *Matter & Energy Cluster child articles, wave 1 · 2026-09-11 · drafted · Physics row P46 · sim pending (matter/Computational_physics).*