Computational Physics

Code Examples

Finite Difference Methods

import numpy as np

def central_difference(f, x, h=1e-5):
    """Central difference for first derivative."""
    return (f(x + h) - f(x - h)) / (2 * h)

def second_derivative(f, x, h=1e-5):
    """Central difference for second derivative."""
    return (f(x + h) - 2*f(x) + f(x - h)) / h**2

def laplacian_2d(f, x, y, h):
    """5-point stencil for 2D Laplacian."""
    return (f(x+h, y) + f(x-h, y) + f(x, y+h) + f(x, y-h) - 4*f(x,y)) / h**2

def solve_poisson_2d(f, nx, ny, Lx, Ly, tol=1e-6):
    """Solve Poisson equation ∇²φ = f using SOR."""
    hx, hy = Lx/nx, Ly/ny
    phi = np.zeros((nx+1, ny+1))
    
    # SOR iteration
    omega = 1.5  # optimal for 2D Poisson
    max_iter = 10000
    
    for iteration in range(max_iter):
        phi_old = phi.copy()
        for i in range(1, nx):
            for j in range(1, ny):
                phi[i,j] = (1-omega)*phi[i,j] + omega/4 * (
                    phi[i+1,j] + phi[i-1,j] + phi[i,j+1] + phi[i,j-1] - hx**2 * f(i*hx, j*hy)
                )
        
        if np.max(np.abs(phi - phi_old)) < tol:
            break
    
    return phi

# Example: Laplace equation on square
def laplace_sor(n=50, tol=1e-6):
    """SOR solver for Laplace equation."""
    phi = np.zeros((n, n))
    phi[:, -1] = 1  # Right boundary = 1
    
    omega = 2 / (1 + np.sin(np.pi/n))
    
    for iteration in range(10000):
        phi_old = phi.copy()
        for i in range(1, n-1):
            for j in range(1, n-1):
                phi[i,j] = (1-omega)*phi[i,j] + omega/4 * (
                    phi[i+1,j] + phi[i-1,j] + phi[i,j+1] + phi[i,j-1]
                )
        
        if np.max(np.abs(phi - phi_old)) < tol:
            print(f"Converged in {iteration} iterations")
            break
    
    return phi

phi = laplace_sor()
print("Laplace equation solved")
print(f"  Max φ = {phi.max():.4f}")

Molecular Dynamics

import numpy as np

class MolecularDynamics:
    def __init__(self, n_atoms, box_size, mass=1.0, dt=0.001):
        self.n = n_atoms
        self.L = box_size
        self.m = mass
        self.dt = dt
        
        # Initialize positions on grid
        n_per_side = int(np.ceil(n_atoms**(1/3)))
        self.r = np.zeros((n_atoms, 3))
        idx = 0
        for i in range(n_per_side):
            for j in range(n_per_side):
                for k in range(n_per_side):
                    if idx < n_atoms:
                        self.r[idx] = [i/n_per_side, j/n_per_side, k/n_per_side] * box_size
                        idx += 1
        
        # Initialize velocities (Maxwell-Boltzmann)
        T = 1.0  # Temperature
        self.v = np.random.normal(0, np.sqrt(T/mass), (n_atoms, 3))
    
    def lennard_jones(self, r):
        """LJ potential: V = 4ε[(σ/r)^12 - (σ/r)^6]"""
        sigma, epsilon = 1.0, 1.0
        r6 = (sigma/r)**6
        return 4 * epsilon * (r6**2 - r6)
    
    def lj_force(self, r_vec):
        """LJ force from potential gradient."""
        sigma, epsilon = 1.0, 1.0
        r = np.linalg.norm(r_vec)
        r2 = r*r
        r6 = (sigma**2/r2)**3
        r12 = r6**2
        force_mag = 24 * epsilon * (2*(sigma**12)/r12 - (sigma**6)/r6) / r2
        return force_mag * r_vec
    
    def compute_forces(self):
        """Calculate forces using neighbor list."""
        f = np.zeros((self.n, 3))
        rc = 2.5  # Cutoff radius
        
        for i in range(self.n):
            for j in range(i+1, self.n):
                r_vec = self.r[i] - self.r[j]
                r_vec = r_vec - self.L * np.round(r_vec / self.L)  # PBC
                r = np.linalg.norm(r_vec)
                
                if r < rc and r > 0:
                    f_ij = self.lj_force(r_vec)
                    f[i] += f_ij
                    f[j] -= f_ij
        
        return f
    
    def integrate(self, n_steps, thermostat=None):
        """Velocity Verlet integration."""
        f = self.compute_forces()
        energies = []
        
        for step in range(n_steps):
            # Update positions
            self.r += self.v * self.dt + 0.5 * f / self.m * self.dt**2
            self.r = self.r % self.L  # Periodic boundaries
            
            # Update forces
            f_new = self.compute_forces()
            
            # Update velocities
            self.v += 0.5 * (f + f_new) / self.m * self.dt
            f = f_new
            
            # Thermostat
            if thermostat:
                thermostat.apply(self.v)
            
            # Energy calculation
            KE = 0.5 * self.m * np.sum(self.v**2)
            PE = self.compute_potential()
            energies.append((KE, PE))
        
        return energies
    
    def compute_potential(self):
        """Calculate total potential energy."""
        PE = 0
        for i in range(self.n):
            for j in range(i+1, self.n):
                r_vec = self.r[i] - self.r[j]
                r_vec = r_vec - self.L * np.round(r_vec / self.L)
                r = np.linalg.norm(r_vec)
                if 0 < r < 2.5:
                    PE += self.lennard_jones(r)
        return PE

md = MolecularDynamics(100, 10.0)
energies = md.integrate(1000)
print(f"Molecular dynamics simulation complete")
print(f"  Initial KE: {energies[0][0]:.2f}")
print(f"  Final KE: {energies[-1][0]:.2f}")

Monte Carlo Methods

import numpy as np
from scipy import integrate

def metropolis_hastings(log_prob, proposal_std, n_samples, x0):
    """Metropolis-Hastings algorithm for sampling."""
    samples = [x0]
    x = x0
    accept_count = 0
    
    for _ in range(n_samples):
        x_proposed = x + np.random.normal(0, proposal_std)
        
        # Acceptance ratio
        log_alpha = log_prob(x_proposed) - log_prob(x)
        
        if np.log(np.random.random()) < log_alpha:
            x = x_proposed
            accept_count += 1
        
        samples.append(x)
    
    return np.array(samples), accept_count / n_samples

# Example: Gaussian mixture sampling
def log_prob_gaussian_mixture(x):
    """Log probability of mixture of two Gaussians."""
    mu1, sigma1 = -2, 1
    mu2, sigma2 = 2, 0.5
    w1, w2 = 0.4, 0.6
    
    p1 = w1 * np.exp(-0.5 * ((x - mu1)/sigma1)**2) / (sigma1 * np.sqrt(2*np.pi))
    p2 = w2 * np.exp(-0.5 * ((x - mu2)/sigma2)**2) / (sigma2 * np.sqrt(2*np.pi))
    return np.log(p1 + p2)

samples, accept_rate = metropolis_hastings(log_prob_gaussian_mixture, 0.5, 10000, 0.0)
print(f"Metropolis-Hastings sampling:")
print(f"  Acceptance rate: {accept_rate*100:.1f}%")
print(f"  Sample mean: {samples.mean():.3f}")
print(f"  Sample std: {samples.std():.3f}")

# Monte Carlo integration
def estimate_pi_mc(n_samples):
    """Estimate π using Monte Carlo integration."""
    np.random.seed(42)
    x = np.random.uniform(-1, 1, n_samples)
    y = np.random.uniform(-1, 1, n_samples)
    inside = np.sum(x**2 + y**2 <= 1)
    return 4 * inside / n_samples

for n in [1000, 10000, 100000, 1000000]:
    pi_est = estimate_pi_mc(n)
    error = abs(pi_est - np.pi)
    print(f"  n = {n:7d}: π ≈ {pi_est:.6f} (error: {error:.6f})")

# Importance sampling
def importance_sampling_integral(f, target_dist, proposal_dist, n_samples):
    """Estimate ∫f(x)p(x)dx using importance sampling."""
    np.random.seed(42)
    x = proposal_dist.rvs(n_samples)
    weights = target_dist.pdf(x) / proposal_dist.pdf(x)
    return np.mean(weights * f(x)), np.std(weights * f(x)) / np.sqrt(n_samples)

from scipy.stats import norm, expon
result, error = importance_sampling_integral(
    lambda x: x**2, 
    norm(0, 1), 
    expon(scale=1), 
    10000
)
print(f"\nImportance sampling estimate: {result:.4f} ± {error:.4f}")

Chaos and Nonlinear Dynamics

import numpy as np
from scipy.integrate import solve_ivp

def lorenz_attractor(t, state, sigma=10, rho=28, beta=8/3):
    """Lorenz system: dx/dt = σ(y-x), dy/dt = x(ρ-z)-y, dz/dt = xy-βz"""
    x, y, z = state
    dxdt = sigma * (y - x)
    dydt = x * (rho - z) - y
    dzdt = x * y - beta * z
    return [dxdt, dydt, dzdt]

def rossler_system(t, state, a=0.2, b=0.2, c=5.7):
    """Rössler system."""
    x, y, z = state
    dxdt = -y - z
    dydt = x + a * y
    dzdt = b + z * (x - c)
    return [dxdt, dydt, dzdt]

def lyapunov_exponent(system, state0, t_span, n_perturbations=4):
    """Estimate Lyapunov exponents using Bennetin's algorithm."""
    t_eval = np.linspace(t_span[0], t_span[1], 100)
    
    # Perturbed states
    perturbed = [state0 + epsilon * np.random.randn(len(state0)) 
                 for epsilon in [1e-8]*n_perturbations]
    
    # Evolve
    for eps_state in perturbed:
        sol = solve_ivp(system, t_span, eps_state, t_eval=t_eval)
    
    return np.random.randn(n_perturbations)  # Placeholder

# Lorenz attractor simulation
state0 = [1.0, 1.0, 1.0]
t_span = (0, 50)
sol = solve_ivp(lorenz_attractor, t_span, state0, t_eval=np.linspace(0, 50, 10000))

print("Lorenz attractor simulation:")
print(f"  Integration complete")
print(f"  x range: [{sol.y[0].min():.2f}, {sol.y[0].max():.2f}]")
print(f"  z range: [{sol.y[2].min():.2f}, {sol.y[2].max():.2f}]")

# Butterfly effect demonstration
state1 = [1.0, 1.0, 1.0]
state2 = [1.0, 1.0, 1.0000001]  # Perturbed

sol1 = solve_ivp(lorenz_attractor, (0, 50), state1, t_eval=np.linspace(0, 50, 10000))
sol2 = solve_ivp(lorenz_attractor, (0, 50), state2, t_eval=np.linspace(0, 50, 10000))

divergence = np.linalg.norm(sol1.y - sol2.y, axis=0)
print(f"\nButterfly effect:")
print(f"  Max separation at t=50: {divergence[-1]:.2f}")
print(f"  Initial separation: 1e-7")

# Period doubling route to chaos
def logistic_map(r, x):
    """x_{n+1} = r x_n (1 - x_n)"""
    return r * x * (1 - x)

print(f"\nLogistic map bifurcation:")
for r in [2.5, 3.0, 3.5, 3.57, 4.0]:
    x = 0.5
    transient = []
    orbit = []
    for _ in range(200):
        x = logistic_map(r, x)
        transient.append(x)
    for _ in range(100):
        x = logistic_map(r, x)
        orbit.append(x)
    
    n_unique = len(set(round(x, 6) for x in orbit))
    behavior = "fixed point" if n_unique == 1 else "period-" + str(n_unique) if n_unique < 10 else "chaotic"
    print(f"  r = {r}: {behavior} ({n_unique} distinct values)")

Finite Element Analysis

import numpy as np

class FiniteElement1D:
    def __init__(self, n_elements, length, degree=1):
        self.n = n_elements
        self.L = length
        self.h = length / n_elements
        
        # Nodes and elements
        self.nodes = np.linspace(0, length, n_elements + 1)
        self.elements = [(i, i+1) for i in range(n_elements)]
    
    def shape_functions(self, xi):
        """Linear shape functions on reference element [-1, 1]."""
        N1 = (1 - xi) / 2
        N2 = (1 + xi) / 2
        return np.array([N1, N2])
    
    def shape_derivatives(self, xi):
        """Derivatives of shape functions."""
        return np.array([-1/2, 1/2])
    
    def assemble_stiffness(self):
        """Assemble stiffness matrix for -u'' = f."""
        K = np.zeros((self.n + 1, self.n + 1))
        
        for (i, j) in self.elements:
            # Element stiffness
            ke = np.array([[1/self.h, -1/self.h], 
                          [-1/self.h, 1/self.h]])
            K[i:i+2, i:i+2] += ke
        
        return K
    
    def assemble_load(self, f):
        """Assemble load vector."""
        F = np.zeros(self.n + 1)
        f_eval = f(self.nodes)
        
        for (i, j) in self.elements:
            fe = self.h / 2 * np.array([f_eval[i], f_eval[j]])
            F[i:i+2] += fe
        
        return F
    
    def solve(self, f, bc_type='dirichlet', bc_values=(0, 1)):
        """Solve boundary value problem."""
        K = self.assemble_stiffness()
        F = self.assemble_load(f)
        
        if bc_type == 'dirichlet':
            # u(0) = 0, u(L) = 1
            K[0, :] = 0
            K[0, 0] = 1
            F[0] = bc_values[0]
            
            K[-1, :] = 0
            K[-1, -1] = 1
            F[-1] = bc_values[1]
        
        return np.linalg.solve(K, F)

# Solve -u'' = 1 with u(0)=0, u(1)=1
def f_constant(x):
    return np.ones_like(x)

fem = FiniteElement1D(100, 1.0)
u = fem.solve(f_constant)

print("FEM solution of -u'' = 1, u(0)=0, u(1)=1:")
print(f"  u(0.5) = {u[50]:.4f} (analytical: 0.375)")
print(f"  u(1.0) = {u[-1]:.4f}")

# Convergence test
def convergence_test():
    errors = []
    ns = [10, 20, 40, 80, 160]
    
    for n in ns:
        fem = FiniteElement1D(n, 1.0)
        u_fem = fem.solve(f_constant)
        u_analytical = 0.5 * fem.nodes - 0.5 * fem.nodes**2
        error = np.max(np.abs(u_fem - u_analytical))
        errors.append(error)
        print(f"  n = {n:3d}: error = {error:.2e}")
    
    # Estimate convergence rate
    if len(errors) > 1:
        rate = np.log(errors[-2]/errors[-1]) / np.log(2)
        print(f"  Convergence rate: ~{rate:.1f}")

convergence_test()

Best Practices

• Verify numerical solutions against known analytical solutions when available to check implementation.
• Use adaptive step sizes in ODE/PDE solvers to balance accuracy and computational cost.
• For molecular dynamics, ensure adequate equilibration before collecting statistics.
• Check energy conservation in symplectic integrators as a validation test.
• For chaotic systems, use many trajectories and ensemble averages for statistical reliability.
• Consider numerical dispersion and dissipation when choosing numerical schemes for wave equations.
• Use appropriate boundary conditions (absorbing, periodic, reflective) for your physical problem.
• Profile and parallelize computationally intensive sections for performance optimization.
• Validate convergence by refining grids/step sizes until results change within tolerance.
• Use dimensional analysis and physical intuition to check if numerical results are reasonable.