Classical Mechanics · Advanced Topics

Computational Physics

Computational physics bridges theory and experiment — solving equations that cannot be solved analytically, simulating systems too complex for closed-form treatment, and testing theoretical predictions with numerical precision. It is now a third pillar of physics alongside theory and experiment.

PrerequisitesDifferentialequations(Ch.DE)⋅Linearalgebra(Ch.LA)⋅Fourieranalysis(Ch.F)⋅LagDifferential equations (Ch. DE) \cdot Linear algebra (Ch. LA) \cdot Fourier analysis (Ch. F) \cdot Lagrangian mechanics (Ch. LA-Mech)
Learning Goals
  • Implement and compare Euler, Verlet, and RK4 integrators and explain why symplectic methods conserve energy over long times.
  • Derive the von Neumann stability condition for explicit finite-difference schemes and state when Crank-Nicolson is preferred.
  • Apply the Metropolis algorithm to sampletheBoltzmanndistributionandidentifytheroleoftheacceptanceprobabilitye(−sample the Boltzmann distribution and identify the role of the acceptance probability e^(-
  • Describe molecular dynamics simulations: force evaluation, Nosé-Hoover thermostats, and the O(N log N) neighbor-list optimisation.
  • Explain spectral (pseudospectral) methods and why they achieve exponential convergence for smooth periodic problems.

CP.1 Numerical Integration of ODEs

Most physics problems reduce to ODEs: ẋ = f(x, t). The simplest integrator:

xn+1=xn+hf(xn,tn)+O(h2)(Euler method, first order)x_{n+1}=x_n+h f(x_n,t_n)+O(h^2) \qquad \text{(Euler method, first order)}(CP.1)

Euler is first-order accurate (local error O(h²), global O(h)). For conservative systems it does not preserve energy. The Leapfrog/Störmer-Verlet method:

xn+1=2xn−xn−1+h2a(xn)(Sto¨rmer-Verlet, symplectic)x_{n+1}=2x_n-x_{n-1}+h^2a(x_n) \qquad \text{(Störmer-Verlet, symplectic)}(CP.2)

Verlet is second-order and symplectic — it preserves the Poincaré invariants (area in phase space) and exhibits bounded energy drift instead of accumulating error. Essential for long-time integrations (molecular dynamics, orbital mechanics).

The Runge-Kutta 4th order (RK4) method is the workhorse for non-conservative problems:

xn+1=xn+h6(k1+2k2+2k3+k4)+O(h5)x_{n+1}=x_n+\frac{h}{6}(k_1+2k_2+2k_3+k_4)+O(h^5)(CP.3)

where k₁ = f(x_n, t_n), k₂ = f(x_n + hk₁/2, t_n + h/2), etc. Fourth-order accurate. Adaptive step-size control: use embedded methods (Runge-Kutta-Fehlberg RK45) that estimate local error and adjust h automatically.

CP.2 Finite Difference Methods for PDEs

Replace continuous derivatives with finite differences on a grid (spacing Δx, Δt):

∂2u∂x2≈ui+1−2ui+ui−1Δx2(second derivative, centered)\frac{\partial^2u}{\partial x^2}\approx\frac{u_{i+1}-2u_i+u_{i-1}}{\Delta x^2} \qquad \text{(second derivative, centered)}(CP.4)

1D Heat equation ∂u/∂t = D ∂²u/∂x²: Explicit FTCS: u_(i,n+1) = uᵢₙ + r(u_(i-1,n) − 2u_(i,n) + u_(i+1,n)) where r = DΔt/Δx². Stability requires r ≤ 1/2 (von Neumann analysis). Implicit Crank-Nicolson: r ≤ ∞ (unconditionally stable, second order in both t and x).

Wave equation ∂²u/∂t² = c² ∂²u/∂x²: FTCS explicit: u_(i,n+1) = 2u_(i,n) − u_(i,n-1) + s²(u_(i+1,n) − 2u_(i,n) + u_(i-1,n)) where s = cΔt/Δx (CFL number). Stability requires s ≤ 1 (Courant-Friedrichs-Lewy condition).

Example CP.1 — Orbital Mechanics with Verlet Integration

Simulate Earth's orbit around the Sun using Störmer-Verlet with adaptive timestep.

Units:Astronomicalunits:1AU=1.496×1011m,1yr=3.156×107s.Intheseunits:GMSun=4π2Astronomical units: 1 AU = 1.496\times10^{11} m, 1 yr = 3.156\times10^{7} s. In these units: G M_{Sun} = 4\pi^{2} AU3/yr2(fromT2=(4π2/GM)a3withT=1yr,a=1AUAU^{3}/yr^{2} (from T^{2} = (4\pi^{2}/GM)a^{3} with T=1yr, a=1AU
Verlet:xn+1=2xn−xn−1+h2a(xn).a=−GM/r3×r(vector).ForEarth:startsat(1,0x_{n+1} = 2x_n - x_{n-1} + h^{2} a(x_{n}). a = -GM/r^{3} \times r (vector). For Earth: starts at (1,0 AUwithvelocity(0,2π)AU/yrAU with velocity (0, 2\pi) AU/yr
Energy conservation:E=12v2−GM/r.WithEuler:Edriftslinearly.WithVerlet:EoscillatesbutremainsboundE = \frac{1}{2}v^{2} - GM/r. With Euler: E drifts linearly. With Verlet: E oscillates but remains bounded—crucialformulti−yearintegrations.ForMercuryprecession:needGRcorrectionaGRed — crucial for multi-year integrations. For Mercury precession: need GR correction a_{GR} = aNewton×(1+3v2/c2+⋯ )toget43′′/centurya_{Newton} \times (1 + 3v^{2}/c^{2} + \cdots) to get 43''/century
Accuracy:Withh=0.01yrandVerlet:perioderror 10−6.Withh=0.01yrandEuler:perioderrorWith h = 0.01 yr and Verlet: period error ~10^{-6}. With h = 0.01 yr and Euler: period error ~10−3.ForthefullSolarSystem(N−body),useREBOUNDorsimilarsymplecticN−bodycode.M10^{-3}. For the full Solar System (N-body), use REBOUND or similar symplectic N-body code. Mcury's orbit unstable on 10^{9} yr timescales (Laskar & Gastineau, 2009.

CP.3 Monte Carlo Methods

Monte Carlo integration: estimate ∫ f(x) dx by sampling random points. For d dimensions, the error scales as N^(−1/2) independent of d — far better than grid methods (which scale as N^(−2/d) in d dimensions, becoming useless for d ≫ 3).

Markov Chain Monte Carlo (MCMC)— Metropolis algorithm: Generate trial move x → x'. Accept if E(x') < E(x). If E(x') > E(x): accept with probability e^(−ΔE/(k_BT)). This samples the Boltzmann distribution P ∝ e^(−E/(k_BT)).

Paccept=min⁡(1,e−ΔE/(kBT))(Metropolis acceptance criterion)P_\mathrm{accept}=\min\left(1,e^{-\Delta E/(k_BT)}\right) \qquad \text{(Metropolis acceptance criterion)}(CP.5)

Applications: the Ising model (compute phase transition), protein conformation sampling, Bayesian posterior sampling, path integrals (lattice QCD computes hadron masses this way).

Quantum Monte Carlo (diffusion Monte Carlo, variational MC): computes exact ground state energies for quantum many-body systems. The Schrödinger equation in imaginary time τ = it: ∂ψ/∂τ = −Hψ → the long-time solution is the ground state (exponentially grows for ground state, decays for excited states). Treat this as a diffusion equation with source/sink from the potential.

CP.4 Molecular Dynamics

Molecular dynamics (MD): integrate Newton's equations for N interacting particles. Force: F_i = −∂U/∂r_i where U = Σ V(rᵢⱼ) (pair potential). Lennard-Jones potential: V(r) = 4ε[(σ/r)¹² − (σ/r)⁶] (repulsion + attraction).

Typical MD: N = 10⁴–10⁷ atoms, timestep h = 1 fs (10⁻¹⁵ s), total time 1 ns–1 μs. The key challenge: force calculation scales as O(N²) naively → reduced to O(N log N) with Verlet neighbor lists and particle-mesh Ewald for electrostatics.

Thermostats control temperature by coupling to a heat bath: Nosé-Hoover (NVT ensemble), Langevin dynamics (adds friction and random force).Barostats control pressure (NPT). Most biological simulations use NPT.

CP.5 Spectral Methods and FFT

For periodic boundary conditions or smooth solutions, expand in Fourier modes. Derivatives become multiplications: (d^n f/dx^n)_k = (ik)^n f_k. The FFT computes all N Fourier coefficients in O(N log N) time.

Pseudospectral method: advance in time in real space (simple); compute derivatives in k-space (accurate). For Navier-Stokes turbulence: dealiasing (2/3 rule removes aliasing error). Spectral accuracy: error ∝ e^(−cN) for smooth functions (exponential convergence vs. algebraic for finite differences).

DMRG (Density Matrix Renormalization Group): O(D³) algorithm to find ground states of 1D quantum systems with matrix product states (MPS). D is the bond dimension. For gapped systems: D grows slowly with system size (area law of entanglement). Extends to 2D (PEPS) and to finite T (MERA, MPO thermofield). State of the art for 1D quantum chemistry.

Definition CP.1 — Common Traps
  • Small timestep is not a proof of correctness: check convergence by reducing h and comparing invariants.
  • Euler is rarely acceptable for long conservative dynamics: energy drift can dominate the physics.
  • Stability and accuracy differ: a stable finite-difference scheme can still be too inaccurate.
  • Random sampling needs equilibration: discard burn-in and account for autocorrelation.
  • Spectral methods assume smoothness: discontinuities cause ringing and destroy exponential convergence.
Exercises — CP.1–CP.5 Computational Physics
1.
Implement one step of RK4 for the nonlinear pendulum \thetä = -(g/L)sin\theta with \theta_{0} = \pi/4, \thetȧ_{0} = 0,h=0.01s.Whatisthesmall−angleperiodT(s)forL=1m0, h = 0.01 s. What is the small-angle period T (s) for L = 1 m
s
Straightforward
2.Write the FTCSfinitedifferenceschemeforthe1Dheatequation.Derivethestabilityconditionr≤FTCS finite difference scheme for the 1D heat equation. Derive the stability condition r \ledoes Crank-Nicolson improve on this?
Intermediate
3.Describe the Metropolis algorithm for the 2D Ising model. What is critical slowing down and how do cluster algorithms solve it?
Intermediate
4.Describeapseudospectralmethodforthe1DBurgersequation\partialu/\partialt+u\partialu/\partialx=ν∂2u/\partialx2.HoDescribe a pseudospectral method for the 1D Burgers equation \partialu/\partialt + u\partialu/\partialx = \nu\partial^{2}u/\partialx^{2}. Howisthestiffviscoustermhandled?Whatistheresolutionrequirementasν→0w is the stiff viscous term handled? What is the resolution requirement as \nu \to 0?
Challenging
Key Takeaways
  • Verlet/leapfrog: symplectic integrator preserving phase space volume. Essential for Hamiltonian systems.
  • RK4: fourth-order, general purpose. Adaptive step-size (RK45) controls global error automatically.
  • FDstability:FTCSheateqrequiresr≤1/2.Waveeq:CFLs≤1.Crank−Nicolson:unconditiFD stability: FTCS heat eq requires r \le 1/2. Wave eq: CFL s \le 1. Crank-Nicolson: unconditionally stable.
  • MetropolisMCMC:acceptalwaysif\DeltaE≤0;accepte−\DeltaE/kTif\DeltaE>0.SamplesBoltzmanMetropolis MCMC: accept always if \DeltaE \le 0; accept e^{-\DeltaE/kT} if \DeltaE > 0. Samples Boltzmann distribution.
  • FFT: O(N log N). Pseudospectral: exponential accuracy for smooth problems + 2/3 dealiasing rule.
  • DMRG:O(D3)groundstatesof1Dquantumsystems.MPS/PEPSaretensornetworkrepresentatioDMRG: O(D^{3}) ground states of 1D quantum systems. MPS/PEPS are tensor network representations.