Computational Physics Scientific¶

What you'll learn¶

• Numerical integration (tested)
• Solving ODEs — simulating motion (tested)
• Monte Carlo methods (tested)
• The accuracy/stability tradeoffs
• The scientific Python stack

Physics is full of equations we can't solve by hand, so we solve them numerically. All three core techniques below are shown in run-verified pure Python.

Numerical integration (tested)¶

Numerical integration in Computational Physics — what it is and when to use it.

Computing an integral (area under a curve) numerically — the trapezoidal rule sums thin trapezoids:

def integrate(f, a, b, n=1000):
    h = (b - a) / n
    total = (f(a) + f(b)) / 2
    for i in range(1, n):
        total += f(a + i * h)
    return total * h

import math
print(round(integrate(lambda x: x*x, 0, 1), 5))    # ∫x² [0,1] = 1/3
print(round(integrate(math.sin, 0, math.pi), 5))    # ∫sin [0,π] = 2

Output:

0.33333
2.0

Both match the exact answers (⅓ and 2). More slices (n) means more accuracy — the essential numerical tradeoff of precision vs computation. This is how you compute work, probability, center of mass — any physical quantity expressed as an integral.

Solving ODEs: simulating motion (tested)¶

Solving ODEs: simulating motion in Computational Physics — what it is and when to use it.

Most physics is differential equations — rates of change. Euler's method steps a system forward in tiny time increments. Here, radioactive decay dy/dt = -y:

def euler(dydt, y0, t0, t1, steps):
    dt = (t1 - t0) / steps
    y, t = y0, t0
    for _ in range(steps):
        y += dydt(t, y) * dt      # step the state forward
        t += dt
    return y

import math
# dy/dt = -y, y(0) = 1  ->  exact answer y(1) = e⁻¹
approx = euler(lambda t, y: -y, y0=1.0, t0=0, t1=1, steps=1000)
print(round(approx, 4), "vs", round(math.exp(-1), 4))

Output:

0.3677 vs 0.3679

Euler's method got 0.3677 against the exact 0.3679 — close, with the small gap being discretization error. Halving the step size roughly halves that error (Euler is first-order). The same technique simulates orbits, pendulums, and circuits — it's the heart of physics engines and the game physics you saw earlier.

Monte Carlo methods (tested)¶

Monte Carlo methods in Computational Physics — what it is and when to use it.

When a problem is too complex to solve directly, Monte Carlo uses random sampling. The classic demo: estimate π by throwing random darts at a square and counting how many land in the inscribed circle:

import random
random.seed(42)

def estimate_pi(n):
    inside = 0
    for _ in range(n):
        x, y = random.random(), random.random()
        if x*x + y*y <= 1:            # inside the quarter circle?
            inside += 1
    return 4 * inside / n             # ratio of areas × 4

print(round(estimate_pi(100_000), 3))

Output:

3.137

100,000 random points estimate π ≈ 3.137 (true value 3.14159). Accuracy improves with more samples — but slowly (error shrinks like 1/√n), which is Monte Carlo's characteristic tradeoff. It shines for high-dimensional integrals and systems (statistical mechanics, quantum Monte Carlo, financial modeling) where deterministic methods break down.

The scientific stack¶

NumPy, SciPy, and Matplotlib as the backbone of numerical physics work.

In real computational physics you'd use optimized libraries:

import numpy as np                    # arrays, vectorized math
from scipy.integrate import solve_ivp, quad   # ODEs and integration

# High-accuracy ODE solve
sol = solve_ivp(lambda t, y: -y, [0, 1], [1.0])

# High-accuracy integration
area, err = quad(lambda x: x**2, 0, 1)   # ~0.3333 with error estimate

Practice exercises¶

1. Increase/decrease n in integrate and observe how the error in ∫x² shrinks.
2. Simulate a falling object with Euler (dv/dt = -g) and compare to the exact v = -gt.
3. Improve the ODE solver to RK4 (midpoint sampling) and compare accuracy to Euler at the same step count.
4. Use Monte Carlo to estimate the area under sin(x) on [0, π] and compare to the tested integrator.
5. Explain why Monte Carlo error scales as 1/√n and what that means for getting one more digit.