Special relativity, quantum and nuclear physics, Monte Carlo, orbit and PDE solvers, fitting with uncertainties, data acquisition and research computing, all in Python.

Computational & Modern Physics

## Why it matters

Many physics problems are too messy for a formula: radiation passing through a detector, neutrons in a reactor, magnets near a phase transition, the path integrals of quantum field theory. **Monte Carlo methods** replace the exact calculation with random sampling. They are easy to write, they work in any number of dimensions, and they come with an honest error bar.

## The core idea

An integral is an average times a size. If `U` is uniform on `[a, b]`, then

```text
integral from a to b of f(x) dx = (b - a) * mean of f(U)
```

So draw `n` random points, average `f` at those points, and multiply by `(b - a)`. The **standard error** of the estimate is

```text
error = (b - a) * s / sqrt(n)
```

where `s` is the sample standard deviation of the `f` values. The key fact: the error shrinks as `1 / sqrt(n)`. Four times as many samples halves the error. That is slow in one dimension, where grid methods do better, but the `1 / sqrt(n)` rate **does not get worse with dimension**. A grid with 10 points per axis in 10 dimensions needs `10^10` evaluations; Monte Carlo does not care.

## Worked example 1: estimating pi

Throw points uniformly in the unit square. The fraction inside the quarter circle `x^2 + y^2 <= 1` approaches `pi / 4`.

```python
import random

rng = random.Random(2026)
n = 100_000
hits = sum(1 for _ in range(n) if rng.random() ** 2 + rng.random() ** 2 <= 1)
print(4 * hits / n)     # about 3.14, typically within 0.01
```

## Worked example 2: an integral with no elementary antiderivative

`integral from 0 to 1 of exp(-x^2) dx = 0.746824...` (it is related to the error function).

```python
import math, random

def mc_integrate(f, a, b, n, seed=0):
    rng = random.Random(seed)
    total = total_sq = 0.0
    for _ in range(n):
        y = f(a + (b - a) * rng.random())
        total += y
        total_sq += y * y
    mean = total / n
    var = (total_sq - n * mean * mean) / (n - 1)
    return (b - a) * mean, (b - a) * math.sqrt(var / n)

est, err = mc_integrate(lambda x: math.exp(-x * x), 0, 1, 100_000, seed=7)
print(f"{est:.4f} +/- {err:.4f}")   # about 0.7468 +/- 0.0006
```

Always report the error with the estimate. If two runs with different seeds disagree by much more than the error bar, something is wrong.

## Worked example 3: a random walk

A particle takes `N` steps of length 1, each in a random direction left or right. Its average position stays 0, but its **mean squared distance** grows linearly: `<x^2> = N`. That is the microscopic picture of diffusion. The typical distance grows only as `sqrt(N)`, which is why a drop of ink spreads quickly at first and slowly later.

```python
rng = random.Random(1)
walkers, steps = 2000, 100
msd = sum(sum(rng.choice((-1, 1)) for _ in range(steps)) ** 2 for _ in range(walkers)) / walkers
print(msd)   # close to 100
```

## Beyond plain sampling

- **Importance sampling:** draw more points where `f` is large to reduce the variance.
- **Markov chain Monte Carlo (Metropolis):** sample from a probability distribution such as the Boltzmann distribution `exp(-E / kT)` by accepting or rejecting random moves. It is the workhorse of statistical physics.
- **Random number quality:** Python's `random` uses the Mersenne Twister, fine for these exercises. NumPy's `np.random.default_rng(seed)` is the modern choice for large simulations.

## Common mistakes

- **Not seeding.** Use a fixed seed for anything you report, and change it on purpose to check stability.
- **Reporting no error bar,** or the standard deviation of `f` instead of the standard error of the mean.
- **Expecting quick precision.** One more correct digit costs 100 times more samples.
- **Reusing the same random numbers** across "independent" runs, which makes results look more certain than they are.

## Recap

- Integral = `(b - a)` times the mean of `f` at uniform random points.
- Error = `(b - a) s / sqrt(n)`; it shrinks as `1 / sqrt(n)` in any dimension.
- Random walks give `<x^2> = N`, the root of diffusion.

## Try it

Open the **Do** tab to write a Monte Carlo integrator that returns an estimate and its standard error.


Estimate integrals and averages by random sampling, quote the statistical error, and model a random walk.

Monte Carlo methods

How to: Monte Carlo Simulation in Python (Introduction)

## Why it matters

Newton's second law turns almost every mechanics problem into a differential equation: `a = F(x) / m`. Only a few cases (a free fall, a spring, a single planet) have neat solutions. For everything else (three bodies, a molecule of a thousand atoms, a galaxy) you step the equations forward in time on a computer. The choice of stepping method decides whether your simulated planet keeps orbiting or spirals away.

## The setup: an orbit in convenient units

Use astronomical units (AU) for distance and years for time. Then the Sun's gravitational parameter is `GM = 4 pi^2 AU^3 / yr^2`, and the Earth's orbit is a circle of radius 1 with speed `2 pi` AU/yr. The acceleration of a planet at position `(x, y)` is

```text
r = sqrt(x^2 + y^2)
a = -GM * (x, y) / r^3
```

## Euler's method and why it fails

The obvious approach updates position and velocity with the current values:

```text
x_new = x + v * dt
v_new = v + a(x) * dt
```

Its error per step is of order `dt^2`, and worse, the errors always push the same way for orbits: energy rises every revolution. With `dt = 0.001` yr for one year, the orbit's energy grows by about 7% and the radius drifts from 1 to about 1.08 AU. Over a thousand years the planet escapes.

## Velocity Verlet

Velocity Verlet uses the acceleration at both ends of the step:

```text
x_new = x + v * dt + 0.5 * a * dt^2
a_new = a(x_new)
v_new = v + 0.5 * (a + a_new) * dt
```

It costs one force evaluation per step, the same as Euler, but its error is second order and it is **symplectic**: it respects the structure of Newtonian mechanics, so the energy wobbles slightly but does not drift. For the same one-year Earth orbit with `dt = 0.001`, the planet returns to within about `1e-4` AU of its start and the relative energy error stays tiny.

```python
import math

GM = 4 * math.pi ** 2

def accel(x, y):
    r3 = (x * x + y * y) ** 1.5
    return -GM * x / r3, -GM * y / r3

x, y, vx, vy, dt = 1.0, 0.0, 0.0, 2 * math.pi, 0.001
ax, ay = accel(x, y)
for _ in range(1000):                 # one year
    x += vx * dt + 0.5 * ax * dt * dt
    y += vy * dt + 0.5 * ay * dt * dt
    nax, nay = accel(x, y)
    vx += 0.5 * (ax + nax) * dt
    vy += 0.5 * (ay + nay) * dt
    ax, ay = nax, nay
print(x, y)                           # close to (1, 0)
```

## Checking a simulation

Never trust an integrator just because it runs. Check quantities that physics says are conserved:

- **Energy per unit mass:** `E = v^2 / 2 - GM / r`. For a circular Earth orbit it is `-2 pi^2`, about `-19.74`.
- **Angular momentum per unit mass:** `L = x * vy - y * vx`, here `2 pi`.
- **Convergence:** halve `dt` and rerun. For a second-order method the error should drop by about 4.

## From two bodies to N bodies

With many bodies, each feels the pull of all the others, so a direct force calculation costs about `N^2` operations per step. Codes add a small **softening length** to avoid infinite forces when two bodies pass very close, and large simulations use tree or mesh methods to approach `N log N`. Molecular dynamics uses exactly the same Verlet step with interatomic forces instead of gravity. Libraries such as SciPy's `solve_ivp` offer adaptive Runge-Kutta methods; they are great for short, accurate runs, but for long conservative simulations a symplectic method like Verlet usually behaves better.

## Choosing the time step

The step must be small compared with the fastest motion in the system. A rule of thumb for orbits is at least a few hundred steps per orbit; for close encounters you need much smaller steps or an adaptive scheme. If halving `dt` changes your answer noticeably, it was too large.

## Common mistakes

- **Using the new position with the old acceleration twice,** which turns Verlet back into a first-order method.
- **Forgetting the `0.5` factors** in the position and velocity updates.
- **Wrong sign on gravity.** Acceleration points toward the Sun: `-GM x / r^3`.
- **Judging by a plot alone.** A drifting energy shows a problem long before the orbit looks wrong.

## Recap

- Euler is simple but drifts; velocity Verlet is just as cheap and conserves energy far better.
- Work in natural units (`AU`, years, `GM = 4 pi^2`) to keep numbers near 1.
- Validate with conserved quantities and a time-step convergence test.

## Try it

Open the **Do** tab to implement velocity Verlet for an orbit and check that energy is conserved.


Integrate a planetary orbit with velocity Verlet, see why Euler's method drifts, and check energy conservation.

Simulating motion

Verlet Integration

## Why it matters

Heat flow, diffusion, waves, electrostatics, fluid flow and the Schrödinger equation are all **partial differential equations** (PDEs): equations in several variables, usually space and time. Most real geometries have no closed-form solution, so physicists discretize them on a grid. The heat equation is the best place to learn this, because it is simple and it shows the main pitfall, numerical instability, very clearly.

## The heat equation

For temperature `u(x, t)` along a thin rod,

```text
du/dt = alpha * d^2u/dx^2
```

`alpha` is the **thermal diffusivity** in m^2/s (about `9.7e-5` for aluminium, `1.1e-4` for copper, about `1.4e-7` for water). The equation says a point warms up when it is colder than the average of its neighbours, so bumps smooth out over time.

## Finite differences

Put the rod on a grid with spacing `dx` and step time by `dt`. Replace each derivative with a difference:

```text
d^2u/dx^2  ~  (u[i+1] - 2*u[i] + u[i-1]) / dx^2
du/dt      ~  (u_new[i] - u[i]) / dt
```

Rearranging gives the **explicit FTCS** update (forward in time, centred in space):

```text
r = alpha * dt / dx^2
u_new[i] = u[i] + r * (u[i+1] - 2*u[i] + u[i-1])
```

With NumPy slicing the whole interior updates in one line:

```python
import numpy as np

def heat_step(u, r):
    new = u.copy()
    new[1:-1] = u[1:-1] + r * (u[2:] - 2 * u[1:-1] + u[:-2])
    return new          # end points unchanged: fixed-temperature boundaries
```

## The stability limit

The explicit scheme is only stable if

```text
r = alpha * dt / dx^2 <= 1/2
```

Break this and small round-off ripples grow exponentially until the numbers explode, even though each step looks harmless. The practical consequence: halving `dx` forces you to cut `dt` by four. For a fine grid, **implicit** methods (backward Euler, Crank-Nicolson) are stable for any `dt` at the cost of solving a linear system each step; SciPy's sparse solvers handle that in real projects.

## Worked example: a rod between two baths

An aluminium rod 0.5 m long starts at 20 °C. At `t = 0` one end is held at 100 °C and the other at 0 °C. With 51 grid points, `dx = 0.01 m`. The largest stable time step is

```text
dt_max = 0.5 * dx^2 / alpha = 0.5 * 1e-4 / 9.7e-5 = 0.52 s
```

Using `dt = 0.4 s` (`r = 0.39`), the profile bends toward a straight line from 100 to 0 °C, which is the **steady state** (where `d^2u/dx^2 = 0`). The slowest-decaying part of the starting error has a time constant of `L^2 / (pi^2 alpha) = 0.25 / (9.87 * 9.7e-5)`, about 260 s, so after about 20 minutes the rod is within roughly a degree of the straight line.

## Testing against an exact solution

Always test a PDE solver on a case you can solve by hand. With both ends held at 0 and a starting profile `u = sin(pi x / L)`, the exact solution keeps its shape and decays:

```text
u(x, t) = sin(pi x / L) * exp(-alpha * pi^2 * t / L^2)
```

Run the solver to some time `t`, compare the peak with `exp(-alpha pi^2 t / L^2)`, and you have a quantitative check. Refining the grid should shrink the difference.

## Common mistakes

- **Ignoring the stability condition** and blaming the physics when values blow up.
- **Updating in place** so new values feed into neighbours in the same step. Compute from a copy (or the slicing form above).
- **Overwriting boundary values** each step. Decide the boundary condition (fixed temperature, insulated, ...) and apply it explicitly.
- **Mixing units,** for example `dx` in cm with `alpha` in m^2/s.
- **Trusting one grid.** Repeat with half the spacing; if the answer moves, refine further.

## Recap

- The heat equation `u_t = alpha u_xx` smooths bumps over time.
- FTCS: `u_new = u + r (u[i+1] - 2u[i] + u[i-1])` with `r = alpha dt / dx^2`.
- Stable only for `r <= 1/2`; implicit methods remove the limit.
- Check solvers against exact solutions and steady states.

## Try it

Open the **Do** tab to build an explicit heat-equation solver with NumPy, including a stability guard and an exact-solution test.


Solve the 1D heat equation with an explicit finite-difference scheme, respect its stability limit and test it against an exact solution.

Solving PDEs

But what is a partial differential equation?  | DE2

Deriving the Heat Equation: A Parabolic Partial Differential Equation for Heat Energy Conservation

Monte Carlo sampling, orbit integration with velocity Verlet and finite-difference solutions of the heat equation.