Physics Simulations in Python

Step 1: Define the Physical System

Most physics problems reduce to ordinary differential equations (ODEs) of the form dy/dt = f(t, y), where y is the state of the system and f describes how the state changes over time. Newton's second law F = ma becomes a system of first-order ODEs by introducing velocity as a separate variable: dx/dt = v, dv/dt = F(x, v, t)/m. In Python, define f as a function: def equations(t, y): x, v = y; F = -k * x; return [v, F / m]. The function takes the current time t and state y, computes forces, and returns the derivatives.

For systems with multiple particles or multiple degrees of freedom, the state vector y contains all positions and velocities. A 2D two-body problem has state [x1, y1, vx1, vy1, x2, y2, vx2, vy2], eight variables. The equations function computes gravitational forces between the bodies and returns all eight derivatives. NumPy arrays make this natural: positions = y[:4].reshape(2, 2), velocities = y[4:].reshape(2, 2), then compute forces using vectorized distance calculations and return the concatenated derivatives.

Conservation laws provide critical validation for simulations. Energy conservation (total kinetic plus potential energy should remain constant in conservative systems), momentum conservation (total momentum is constant in isolated systems), and angular momentum conservation each give you a quantity that you can compute at every time step and verify remains constant. If the conserved quantity drifts significantly during the simulation, the numerical method is accumulating errors and you need a smaller time step, a higher-order solver, or a symplectic integrator that preserves conservation laws by construction.

Step 2: Solve with Numerical Integration

scipy.integrate.solve_ivp is the standard solver for initial value ODE problems. from scipy.integrate import solve_ivp. sol = solve_ivp(equations, [0, 10], y0, t_eval=np.linspace(0, 10, 1000), method='RK45', rtol=1e-8, atol=1e-10). t_eval specifies the output times. rtol and atol control the relative and absolute error tolerances: tighter tolerances produce more accurate results at the cost of more computation. The default RK45 (4th/5th order Runge-Kutta) handles most non-stiff problems well.

Solver choice depends on the problem characteristics. RK45 is the default for smooth, non-stiff problems (projectile motion, orbital mechanics, simple oscillators). RK23 is less accurate but faster for rough solutions. DOP853 is an 8th-order method for problems requiring high accuracy (precision orbital mechanics, long-duration simulations). BDF and Radau handle stiff problems where some timescales are much shorter than others (chemical kinetics with fast and slow reactions, electrical circuits with widely separated time constants). If RK45 takes extremely long or produces warnings about excessive steps, switch to BDF or Radau.

For long-duration simulations where energy conservation matters (planetary orbits over thousands of years, molecular dynamics), symplectic integrators preserve the geometric structure of Hamiltonian systems. The Verlet method (also called leapfrog or Stormer-Verlet) is the simplest: x_new = 2*x - x_old + a * dt**2 for position, v = (x_new - x_old) / (2 * dt) for velocity. Symplectic integrators do not minimize error at each step (they are only second-order accurate), but they prevent the long-term energy drift that causes non-symplectic methods to produce unphysical results over many orbits. Implement Verlet in a simple loop: it is only a few lines of code and produces qualitatively correct long-term behavior that higher-order non-symplectic methods cannot match.

Step 3: Simulate Classical Mechanics

Projectile motion with air resistance demonstrates the basic workflow. Without air resistance, the trajectory is a parabola with an analytical solution. With quadratic air drag (F_drag = -0.5 * rho * C_d * A * v * |v|), no closed-form solution exists and numerical simulation is required. Define: def projectile(t, y): x, z, vx, vz = y; v = np.sqrt(vx**2 + vz**2); drag = 0.5 * rho * Cd * A * v; return [vx, vz, -drag * vx / m, -g - drag * vz / m]. Solve and plot x vs z to see how air resistance shortens the range and lowers the peak height compared to the ideal parabola.

Coupled oscillators model everything from molecules to bridges. Two masses connected by springs: m1*a1 = -k1*x1 + k_c*(x2 - x1), m2*a2 = -k2*x2 + k_c*(x1 - x2), where k_c is the coupling spring constant. The system exhibits normal modes: collective oscillation patterns where all masses move at the same frequency. For N coupled oscillators, the eigenvalues of the stiffness matrix give the normal mode frequencies, and the eigenvectors give the mode shapes. Simulate and visualize the energy transferring between oscillators when the coupling is weak, or the coherent motion of normal modes when you initialize in an eigenmode.

N-body gravitational simulation models planetary systems, star clusters, and galaxy dynamics. Each body exerts gravitational force on every other body: F_ij = G * m_i * m_j / r_ij**2 in the direction from i to j. The total force on each body is the sum over all other bodies. For N bodies, this is O(N^2) force calculations per time step. Simulate the inner solar system (Sun plus 4 planets) by using real masses, distances (in AU), and velocities (from NASA JPL Horizons). Verify that your simulation reproduces known orbital periods. For larger N (star clusters with thousands of bodies), use the Barnes-Hut algorithm (O(N log N)) or fast multipole methods.

Chaotic systems demonstrate sensitive dependence on initial conditions. The double pendulum is the classic example: two connected pendulum arms whose motion is deterministic but chaotic, meaning tiny differences in initial angle grow exponentially over time until the trajectories diverge completely. Simulate two double pendulums with initial angles differing by 0.001 radians and plot both trajectories: they track each other initially, then diverge into completely different patterns. Compute the Lyapunov exponent (the rate of divergence) to quantify the chaos. The Lorenz system (three coupled ODEs) is another classic chaotic system whose strange attractor produces the famous butterfly-shaped trajectory.

Step 4: Simulate Fields and Waves

Partial differential equations (PDEs) require discretizing space as well as time. The finite difference method replaces spatial derivatives with differences on a grid. For the 1D heat equation (du/dt = alpha * d2u/dx2): create a spatial grid x = np.linspace(0, L, N), approximate the second derivative as (u[i+1] - 2*u[i] + u[i-1]) / dx**2, and step forward in time with the explicit Euler method: u_new[i] = u[i] + alpha * dt / dx**2 * (u[i+1] - 2*u[i] + u[i-1]). The stability condition dt < dx**2 / (2 * alpha) must be satisfied or the simulation explodes. Plot temperature profiles at different times to see heat diffusing from hot regions to cold.

Wave equation simulation (d2u/dt2 = c**2 * d2u/dx2) uses the same spatial discretization but requires a second-order time stepping: u_new[i] = 2*u[i] - u_old[i] + (c * dt / dx)**2 * (u[i+1] - 2*u[i] + u[i-1]). The Courant number (c * dt / dx) must be less than or equal to 1 for stability. Initialize with a Gaussian pulse and watch it split into two pulses traveling in opposite directions. Add boundary conditions: fixed boundaries (u = 0) cause reflections with phase inversion, free boundaries (du/dx = 0) cause reflections without inversion, and absorbing boundaries minimize reflections at domain edges.

Two-dimensional simulations use 2D grids and five-point stencils for the Laplacian: nabla2_u[i,j] = (u[i+1,j] + u[i-1,j] + u[i,j+1] + u[i,j-1] - 4*u[i,j]) / dx**2. The 2D heat equation produces diffusion patterns. The 2D wave equation produces circular wave fronts from point sources, interference patterns from multiple sources, and diffraction around obstacles. Visualize with plt.imshow(u, cmap='hot') for temperature fields or plt.imshow(u, cmap='RdBu', vmin=-1, vmax=1) for wave amplitude fields. Animate with matplotlib.animation.FuncAnimation to show the temporal evolution.

Electrostatic potential problems solve Laplace's equation (nabla2_V = 0) or Poisson's equation (nabla2_V = -rho/epsilon_0) with boundary conditions specifying voltages on conductors. The iterative Jacobi or Gauss-Seidel method repeatedly updates each grid point to the average of its neighbors until convergence. For a parallel plate capacitor: set V = +1 on one plate, V = -1 on the other, iterate until the maximum change per iteration drops below a threshold, then compute the electric field as E = -gradient(V). Plot field lines with plt.streamplot() and equipotential contours with plt.contour().

Step 5: Visualize and Animate Results

Phase space plots reveal the qualitative behavior of dynamical systems. Plot velocity vs position (v vs x) for oscillators: simple harmonic motion produces ellipses, damped motion produces inward spirals, and driven motion at resonance produces expanding spirals that stabilize into limit cycles. For the double pendulum, plot angle1 vs angle2 to see the chaotic attractor. Phase portraits show the topology of the dynamics: fixed points (equilibria), limit cycles (periodic orbits), and strange attractors (chaotic behavior) are all visible as geometric structures in phase space.

Energy plots verify simulation accuracy and reveal physical behavior. Plot kinetic energy, potential energy, and total energy vs time on the same axes. For conservative systems, total energy should remain constant (any drift indicates numerical error). For damped systems, total energy should decrease monotonically. For driven systems, energy oscillates and may grow (resonance) or stabilize (steady state). The energy plot is the first diagnostic you should check after every simulation: if energy conservation fails, the simulation results cannot be trusted.

Animations show temporal evolution that static plots cannot capture. from matplotlib.animation import FuncAnimation. Create a figure and an initial plot, then define an update function that modifies the plot data for each frame. ani = FuncAnimation(fig, update, frames=range(num_frames), interval=30, blit=True). Save as video: ani.save('simulation.mp4', writer='ffmpeg', dpi=150). Save as GIF: ani.save('simulation.gif', writer='pillow', fps=30). Animations are essential for wave simulations, fluid dynamics, and any system where the spatial pattern evolves in time.

3D visualization with matplotlib's Axes3D or with Mayavi shows trajectories in 3D space (planetary orbits, particle tracks), surfaces (potential energy surfaces, wave functions), and vector fields (electromagnetic fields, fluid velocity). ax = fig.add_subplot(projection='3d'). ax.plot3D(x, y, z) for trajectories. ax.plot_surface(X, Y, Z, cmap='viridis') for surfaces. For publication, 3D plots should include clear axis labels, reasonable viewing angles that reveal the important structure, and possibly multiple views (top, side, perspective) in a multi-panel figure.