A mid-level course for future physicists that joins modern physics to the numerical methods physicists use every day. You will compute time dilation and the Lorentz factor without losing precision; find the energy levels and wavefunctions of a particle in a box; model radioactive decay both exactly and by simulation; estimate integrals with Monte Carlo sampling and quote their error; integrate an orbit with the velocity Verlet method and check energy conservation; solve the heat equation with finite differences and respect the stability limit; fit data with uncertainties and judge the fit with chi-squared; record and analyse real sensor data with a phone or a free microcontroller simulator; and see how clusters, schedulers and reproducible research software work. Prerequisites: Scientific Python and Physics for Engineering & Science. Code activities run Python in your browser with NumPy where it helps. This course teaches knowledge and skills; physicist roles need at least a bachelor's degree in physics, and most research positions in national labs, universities and R&D need a PhD.

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

Special relativity is not only for astrophysics. Particle accelerators, cosmic-ray detectors and GPS receivers all need it to give correct answers. Muons made in the upper atmosphere reach the ground only because their clocks run slow in our frame, and GPS satellite clocks would drift by microseconds a day without relativistic corrections. As a physicist you will compute these effects often, and the calculation has a numerical trap that is worth learning early.

## The core idea

Einstein's two postulates are simple: the laws of physics are the same in every inertial (non-accelerating) frame, and light in a vacuum travels at the same speed `c = 299,792,458 m/s` for every observer. Everything else follows from taking the second one seriously.

The key quantity is the **Lorentz factor**:

```text
beta  = v / c
gamma = 1 / sqrt(1 - beta^2)
```

- `gamma` is 1 at rest and grows without limit as `v` approaches `c`.
- **Time dilation:** a clock moving at speed `v` ticks slowly. If a process takes the *proper time* `t0` in its own rest frame, we measure `t = gamma * t0`.
- **Length contraction:** an object of rest length `L0` is measured as `L = L0 / gamma` along the direction of motion.
- **Energy:** total energy is `E = gamma * m * c^2`, and kinetic energy is `(gamma - 1) * m * c^2`.

At 60% of light speed, `beta = 0.6`, `1 - 0.36 = 0.64`, `sqrt(0.64) = 0.8`, so `gamma = 1.25`. At 80%, `gamma = 1 / 0.6 = 5/3`.

## Worked example: cosmic-ray muons

A muon has a mean lifetime of 2.197 µs in its own rest frame. Muons are created about 10 to 15 km up and many travel at around `0.995 c`.

```text
beta^2 = 0.995^2 = 0.990025
gamma  = 1 / sqrt(0.009975) = 1 / 0.099875 = 10.01
t (lab) = 10.01 * 2.197 µs = 22.0 µs
distance in one lifetime = 0.995 * 2.998e8 m/s * 22.0e-6 s = 6.56 km
```

Without time dilation they would travel only about 655 m in one lifetime, so almost none would reach a detector at sea level. In the muon's own frame the explanation is length contraction instead: the atmosphere is about ten times thinner. Both frames agree on what the detector counts.

## Worked example: the round-off trap

GPS satellites move at about 3.9 km/s. Here `beta^2` is about `1.7e-10`, so `gamma - 1` is about `8.5e-11`. A double-precision float carries about 16 significant digits, and `1 - beta^2` sits so close to 1 that most of those digits are spent on the leading 1. Subtracting 1 from `gamma` afterwards throws away the digits you wanted.

```python
import math

c = 299_792_458.0
v = 300.0                           # an airliner, in m/s
b2 = (v / c) ** 2                   # about 1.0e-12

naive = 1 / math.sqrt(1 - b2) - 1   # only a few correct digits
stable = math.expm1(-0.5 * math.log1p(-b2))  # full precision
series = b2 / 2 + 3 * b2**2 / 8     # Taylor series check
print(naive, stable, series)
```

`math.log1p(x)` computes `ln(1 + x)` and `math.expm1(x)` computes `e^x - 1`, both accurately for tiny `x`. Since `gamma = (1 - b2)^(-1/2) = exp(-0.5 * ln(1 - b2))`, the stable form gives `gamma - 1` directly. This pattern (rewrite to avoid subtracting nearly equal numbers) appears all over computational physics.

For scale: the satellite's speed makes its clock lose about 7 µs per day, while the weaker gravity at orbit (a general relativity effect) makes it gain about 45 µs per day. The net is about +38 µs per day, which would cause kilometres of position error a day if not corrected.

## Common mistakes

- **Mixing up which time is proper time.** The proper time is measured by a clock at rest relative to the process (the muon's own clock). The dilated time is always the longer one.
- **Forgetting units.** Keep `v` and `c` in the same units, or work in `beta` from the start.
- **Allowing `v >= c`.** A massive object can't reach `c`; your code should reject it rather than return `inf` or a complex number.
- **Using `gamma - 1` computed naively** for slow speeds, then trusting the digits.
- **Expecting contraction across the motion.** Only lengths along the direction of motion contract.

## Recap

- `gamma = 1 / sqrt(1 - v^2/c^2)`; `t = gamma * t0`; `L = L0 / gamma`; `E = gamma * m * c^2`.
- Muons reaching the ground are direct evidence of time dilation.
- For small speeds, compute `gamma - 1` with `log1p`/`expm1` or a series, never by subtracting.

## Try it

Open the **Do** tab and write a small relativity toolkit, including a precise `gamma - 1` for slow speeds.


Compute the Lorentz factor, time dilation and length contraction, and avoid the round-off trap at everyday speeds.

Special relativity

Special Relativity Part 2: Time Dilation and the Twin Paradox

## Why it matters

Quantum mechanics explains why atoms are stable, why materials conduct or insulate, why LEDs and quantum dots have the colours they do, and how lasers, transistors and MRI work. You do not need the full theory to start computing with it. The **particle in a box** is the standard first model: simple enough to solve by hand and rich enough to show every key idea.

## The core ideas

- **Wavefunction.** A particle is described by a wavefunction `psi(x)`. It is not a physical wave you can see; it is a tool for predictions.
- **Probability.** `|psi(x)|^2` is a probability density. The chance of finding the particle between `a` and `b` is the integral of `|psi|^2` from `a` to `b`. The total over all space must be 1 (the wavefunction is **normalized**).
- **Schrödinger equation.** For a time-independent problem in one dimension:

```text
-(hbar^2 / 2m) * psi''(x) + V(x) * psi(x) = E * psi(x)
```

Only some values of `E` give wavefunctions that obey the boundary conditions. That is why energy is **quantized**.

## The particle in a box

Put a particle of mass `m` in a region `0 < x < L` with zero potential inside and infinitely high walls outside. The wavefunction must be zero at both walls. Sine waves that fit a whole number of half-wavelengths into the box are the solutions:

```text
psi_n(x) = sqrt(2/L) * sin(n * pi * x / L),   n = 1, 2, 3, ...
E_n = n^2 * h^2 / (8 * m * L^2)
```

What this tells you:

- Energy grows as `n^2`: levels 1, 4, 9, 16 times the ground state.
- The lowest energy is not zero. Confining a particle always costs energy (**zero-point energy**), a consequence of the uncertainty principle.
- State `n` has `n - 1` nodes inside the box where the particle is never found.
- A smaller box or a lighter particle spreads the levels further apart (`1/(m L^2)`).

## Worked example: an electron in a 1 nm box

```text
h   = 6.62607015e-34 J s
m_e = 9.1093837e-31 kg
L   = 1.0e-9 m

E_1 = h^2 / (8 m L^2) = 4.390e-67 / 7.288e-48 = 6.025e-20 J
    = 6.025e-20 / 1.602e-19 = 0.376 eV
```

So `E_2 = 4 * 0.376 = 1.504 eV` and `E_3 = 3.384 eV`. If the electron drops from `n = 3` to `n = 2`, it releases `5 * 0.376 = 1.880 eV`, and the photon wavelength is

```text
lambda = h c / dE = 1239.84 eV nm / 1.880 eV = 659 nm   (red light)
```

This is roughly how colour works in dye molecules and quantum dots: the size of the region the electron is confined to sets the colour. Real atoms have a different potential, so their levels follow other formulas. For hydrogen, `E_n = -13.6 eV / n^2`, and the 3 to 2 transition gives 656 nm, the red H-alpha line.

## Worked example: probabilities

For the ground state, what is the chance of finding the particle in the middle half of the box?

```text
P = integral from L/4 to 3L/4 of (2/L) sin^2(pi x / L) dx = 1/2 + 1/pi = 0.818
```

A classical particle bouncing at constant speed would spend exactly 50% of its time there. The quantum ground state piles probability up in the middle. Many such integrals have no neat closed form, so you compute them numerically:

```python
import math

def psi(n, x, L):
    return math.sqrt(2 / L) * math.sin(n * math.pi * x / L)

def prob(n, L, a, b, steps=10_000):
    dx = (b - a) / steps
    return sum(psi(n, a + (i + 0.5) * dx, L) ** 2 for i in range(steps)) * dx

print(prob(1, 1.0, 0.25, 0.75))   # 0.8183...
```

This is the midpoint rule: split the interval, evaluate at each centre, add up.

## Common mistakes

- **Treating `psi` as the probability.** The probability density is `|psi|^2`, and probabilities come from integrating it.
- **Starting `n` at 0.** `n = 0` gives `psi = 0` everywhere, which is no particle at all.
- **Mixing `h` and `hbar`.** The box formula above uses `h`. With `hbar = h / (2 pi)` it reads `n^2 pi^2 hbar^2 / (2 m L^2)`; both give the same number.
- **Forgetting the eV conversion** (`1 eV = 1.602176634e-19 J`) and quoting joules as if they were electronvolts.

## Recap

- `|psi|^2` is a probability density; normalized wavefunctions integrate to 1.
- Boundary conditions quantize energy: `E_n = n^2 h^2 / (8 m L^2)` in a box.
- Transitions emit photons with `lambda = h c / dE`.
- Numerical integration (midpoint rule) handles the probabilities you can't do by hand.

## Try it

Open the **Do** tab to compute box energy levels, a transition wavelength and a probability by numerical integration.


Explain wavefunctions, probability and quantized energy levels using the particle in a box, and compute them in Python.

Quantum basics

Quantum Mechanics - Part 1: Crash Course Physics #43

Particle in a Box Part 1: Solving the Schrödinger Equation

## Why it matters

Nuclear physics underlies carbon dating, medical imaging and cancer therapy, smoke detectors, nuclear power and the detectors used in particle and astro physics. Radioactive decay is also the cleanest example of a **random process with a precise average**, which makes it the perfect bridge from physics to the simulation methods in the next module.

## Atoms and nuclei

An atom is a tiny, dense nucleus surrounded by electrons. The nucleus holds `Z` protons and `N` neutrons; the mass number is `A = Z + N`. Isotopes of an element share `Z` but differ in `N`, for example carbon-12 (6 protons, 6 neutrons) and carbon-14 (6 protons, 8 neutrons).

The strong nuclear force holds the nucleus together against the electric repulsion between protons. A nucleus weighs **less** than its separate protons and neutrons. The missing mass is the **binding energy**, through `E = m c^2`, with `1 u = 931.494 MeV/c^2`.

**Worked example: the deuteron.** A deuteron is one proton plus one neutron.

```text
proton  = 1.007276 u
neutron = 1.008665 u
sum     = 2.015941 u
deuteron= 2.013553 u
mass defect = 0.002388 u  ->  0.002388 * 931.494 = 2.224 MeV
```

That is millions of times the few eV involved in chemical bonds, which is why nuclear reactions release so much energy.

## Kinds of radiation

- **Alpha:** a helium-4 nucleus (2 protons, 2 neutrons). Stopped by a sheet of paper or the outer layer of skin, but dangerous if the source is inhaled or swallowed.
- **Beta:** an electron or positron from a neutron turning into a proton (or the reverse). Stopped by a few mm of aluminium or plastic.
- **Gamma:** a high-energy photon. Needs dense shielding such as lead or concrete.

## The decay law

Each unstable nucleus has the same fixed chance of decaying per unit time, whatever its age. That single fact gives exponential decay:

```text
N(t) = N0 * exp(-lambda * t)
lambda = ln(2) / T_half
activity A = lambda * N       (decays per second, becquerel: 1 Bq = 1 decay/s)
```

After one half-life half remains, after two a quarter, after `k` half-lives `(1/2)^k`. The older unit, the curie, is `1 Ci = 3.7e10 Bq`.

**Worked example: carbon dating.** Living things keep a steady ratio of carbon-14 to carbon-12. After death the C-14 decays with `T_half` of about 5,730 years (the conventional value; recent measurements give about 5,700). A sample with 25% of the living ratio has been through two half-lives: about 11,460 years.

**Worked example: activity.** Iodine-131 (used in thyroid treatment) has `T_half = 8.02 days`.

```text
lambda = 0.693147 / (8.02 * 86,400 s) = 1.0003e-6 per second
A = 1.0003e-6 * 1.0e15 atoms = 1.0e9 Bq (1 GBq)
```

## Simulating decay

Analytic formulas give the average. Experiments see **counts that fluctuate** around it. To model that, give every nucleus a decay probability `p = 1 - exp(-lambda * dt)` per time step and roll a random number for each one:

```python
import math, random

def simulate(n0, half_life, dt, steps, seed=1):
    rng = random.Random(seed)            # fixed seed = reproducible run
    p = 1 - math.exp(-math.log(2) / half_life * dt)
    n, counts = n0, [n0]
    for _ in range(steps):
        n -= sum(1 for _ in range(n) if rng.random() < p)
        counts.append(n)
    return counts
```

With 20,000 nuclei the result after one half-life will be close to 10,000 but not exactly. The spread follows counting statistics: a count of about `N` has a standard deviation near `sqrt(N)`, so about ±100 here. A good simulation agrees with theory within that noise, not to the last digit.

## Safety and real work

Radioactive sources, X-ray equipment and accelerators are regulated. In the US, licences come from the Nuclear Regulatory Commission or a state agency, and every lab has a radiation safety officer and required training. Never collect or handle radioactive materials at home; everything in this course is done with simulations and published data.

## Common mistakes

- **Using `lambda = 1 / T_half`.** The decay constant is `ln(2) / T_half`.
- **Mixing time units** (half-life in days, `t` in seconds). Convert first.
- **Expecting a simulation to match the formula exactly.** Compare within the expected statistical spread.
- **Leaving the random seed unset** in results you report, so nobody (including you) can reproduce them.

## Recap

- `A = Z + N`; binding energy is mass defect times `c^2`.
- `N(t) = N0 exp(-lambda t)`, `lambda = ln 2 / T_half`, `A = lambda N`.
- Simulated decays fluctuate by about `sqrt(N)`; seed your generator.

## Try it

Open the **Do** tab to code the decay law and a seeded Monte Carlo decay simulation, then compare them.


Describe nuclei, binding energy and radiation, and model radioactive decay both exactly and by simulation.

Atomic and nuclear physics

Nuclear Half Life: Intro and Explanation

Half life | Radioactivity | Physics | FuseSchool

Special relativity, the quantum particle in a box and atomic and nuclear physics, each turned into a small Python model.

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.

Computational Methods

## Why it matters

A measurement without an uncertainty is not a result. Physicists fit models to data every day: a decay constant from counts, a spring constant from extensions, `g` from a pendulum. The fit gives the best parameters, but the questions that matter are **how sure are we** and **does the model actually fit**. Weighted least squares and chi-squared answer both.

## Uncertainties first

Each data point `y_i` has an uncertainty `sigma_i`, its standard deviation. Points with small `sigma` should count more. Two rules for propagating uncertainties you will use constantly (for independent errors):

```text
sum or difference:  sigma_f^2 = sigma_a^2 + sigma_b^2
product or ratio:   (sigma_f / f)^2 = (sigma_a / a)^2 + (sigma_b / b)^2
```

## Weighted least squares for a line

To fit `y = m x + b`, choose `m` and `b` that minimise

```text
chi^2 = sum of ((y_i - (m x_i + b)) / sigma_i)^2
```

With weights `w_i = 1 / sigma_i^2`, the answer has a closed form:

```text
S = sum w     Sx = sum w x     Sy = sum w y
Sxx = sum w x^2     Sxy = sum w x y     D = S*Sxx - Sx^2

m = (S*Sxy - Sx*Sy) / D          sigma_m = sqrt(S / D)
b = (Sxx*Sy - Sx*Sxy) / D        sigma_b = sqrt(Sxx / D)
```

## Judging the fit: chi-squared

Each term in `chi^2` asks: how many sigmas away is this point from the line? If the model is right and the `sigma` values are honest, each term is about 1 on average. With `N` points and 2 fitted parameters there are `N - 2` **degrees of freedom**, and

```text
reduced chi^2 = chi^2 / (N - 2)   ~ 1 for a good fit
```

- **Much larger than 1:** the model is wrong, or the uncertainties are underestimated.
- **Much smaller than 1:** the uncertainties are probably overestimated.
- Always look at the **residuals** too: a pattern (a curve, a trend) means the model misses something even if `chi^2` looks fine.

## Worked example: a cart on a track

Position of a cart every second, each with `sigma = 0.03 m`:

```text
t (s): 0     1     2     3     4     5
x (m): 0.02  0.51  1.03  1.48  2.05  2.49
```

The weighted fit gives

```text
m = 0.498 +/- 0.007 m/s     (the speed)
b = 0.019 +/- 0.022 m       (consistent with starting at 0)
chi^2 = 3.60 for 4 degrees of freedom -> reduced chi^2 = 0.90
```

The reduced chi-squared near 1 says a constant speed describes the data at this precision. Report it as `v = 0.498 +/- 0.007 m/s`: round the uncertainty to one or two significant figures and the value to match.

## With NumPy and SciPy

```python
import numpy as np

t = np.array([0, 1, 2, 3, 4, 5.0])
x = np.array([0.02, 0.51, 1.03, 1.48, 2.05, 2.49])
sigma = np.full(6, 0.03)

coef, cov = np.polyfit(t, x, 1, w=1 / sigma, cov="unscaled")
print(coef, np.sqrt(np.diag(cov)))   # [0.498 0.019] [0.0072 0.0217]
```

Note that `np.polyfit` takes weights `1/sigma`, not `1/sigma^2`, and `cov="unscaled"` returns errors based on your sigmas. For nonlinear models, `scipy.optimize.curve_fit(f, x, y, sigma=sigma, absolute_sigma=True)` does the same job; `absolute_sigma=True` tells it your sigmas are real uncertainties rather than relative weights.

## Common mistakes

- **Unweighted fits for data with different uncertainties,** which lets noisy points pull the line.
- **Passing `1/sigma^2` to `polyfit`** (it squares the weights itself).
- **Quoting too many digits,** such as `0.4977143 +/- 0.0071714`.
- **Ignoring a reduced chi-squared of 5** because the line "looks fine".
- **Forgetting the degrees of freedom:** divide by `N - p`, where `p` is the number of fitted parameters.

## Recap

- Weighted least squares minimises chi-squared with weights `1/sigma^2`.
- Parameter errors come from the same sums: `sigma_m = sqrt(S/D)`, `sigma_b = sqrt(Sxx/D)`.
- Reduced chi-squared near 1 plus patternless residuals means a believable fit.

## Try it

Open the **Do** tab to write a weighted line fit with parameter uncertainties and chi-squared, checked against NumPy.


Fit a straight line to data with uncertainties, quote the parameter errors, and judge the fit with chi-squared.

Data analysis and fitting

Least Squares Fit and Graphing in Python

Chi Squared

Install phyphox (or open a new Arduino Uno project in Wokwi and add a potentiometer wired to 5V, GND and A0, with the logger sketch).

Lay the phone flat and still on a table. In phyphox open 'Acceleration with g' and record for 30 seconds without touching it (in Wokwi, run the sketch for 30 seconds without moving the potentiometer).

Export the data as CSV (phyphox: menu, Export Data, CSV) or copy the serial output into a text file. Write down the date, device model and settings in a notes file next to it.

Load the file in Python with pandas or numpy. Find the sample rate from the timestamps (number of samples divided by elapsed time) and check whether the spacing is steady.

Compute the mean and standard deviation of the z-axis (or voltage) column. The phone's mean should be close to 9.81 m/s^2; note how far off it is.

Compute the standard error of the mean (std / sqrt(N)) and compare it with the difference from 9.81. Decide whether the difference is noise or a calibration offset.

Repeat with the phone resting on its long edge (or turn the potentiometer to a new fixed position) and check that the reading moves to a different axis (or voltage) as expected.

Write three sentences: the sample rate, the noise level per sample, and one systematic effect you would calibrate out in a real experiment.

## Why it matters

Every physics result starts as a voltage in an instrument. **Data acquisition (DAQ)** is the chain that turns a physical quantity (temperature, acceleration, light, charge) into numbers in a file. Physicists design, calibrate and debug these chains constantly, and many analysis problems are really acquisition problems: the wrong sample rate, a missing timestamp, an uncalibrated sensor.

## The DAQ chain

1. **Sensor or transducer:** converts the quantity into an electrical signal. Examples: a thermistor (resistance changes with temperature), a photodiode (current from light), a MEMS accelerometer (capacitance changes with acceleration), a photomultiplier (charge pulses from single photons).
2. **Signal conditioning:** amplifiers, filters and dividers that bring the signal into the converter's range and remove noise. An **anti-aliasing low-pass filter** removes frequencies the converter can't sample correctly.
3. **Analog-to-digital converter (ADC):** samples the voltage at fixed times and rounds it to an integer.
4. **Timing and storage:** a clock stamps every sample; software writes the data with its metadata.

## Resolution

An `n`-bit ADC splits its input range into `2^n` steps.

```text
step size = V_range / 2^n
10-bit, 0-5 V:   5 / 1024  = 4.9 mV per count
12-bit, 0-3.3 V: 3.3 / 4096 = 0.81 mV per count
16-bit, +/-10 V: 20 / 65536 = 0.31 mV per count
```

A reading of 512 counts on a 10-bit, 5 V ADC is `512 * 5 / 1024 = 2.50 V`. Resolution is not accuracy: an ADC with tiny steps can still be offset or noisy, which is why you **calibrate** against a known reference.

## Sampling rate and aliasing

The **Nyquist criterion** says you must sample faster than twice the highest frequency in the signal. A 60 Hz signal sampled at 100 samples per second does not just look rough, it masquerades as a 40 Hz signal (`100 - 60`). Once aliased, no software can undo it, so you filter before the ADC and choose a sample rate well above twice the frequencies you care about (five to ten times is common).

## Noise and averaging

Random noise shrinks when you average: the standard error of `N` independent readings is `sigma / sqrt(N)`. Averaging does nothing for systematic errors, such as a sensor that reads 0.2 °C high; only calibration fixes those.

## Phones and microcontrollers as DAQ

- **Smartphone sensors** are real MEMS instruments. The free **phyphox** app, from RWTH Aachen University (https://phyphox.org), records the accelerometer, gyroscope, magnetometer, light sensor and microphone and exports CSV files.
- **Microcontrollers** such as Arduino boards read sensors through built-in ADCs (`analogRead` returns 0 to 1023 on a 10-bit, 5 V Uno) and send data over USB serial. The free **Wokwi** browser simulator (https://wokwi.com) runs Arduino code with virtual sensors and potentiometers, so you can practise without hardware.

A minimal Arduino logger looks like this:

```cpp
void setup() { Serial.begin(115200); }
void loop() {
  int counts = analogRead(A0);             // 0..1023
  float volts = counts * 5.0 / 1023.0;
  Serial.print(millis()); Serial.print(",");
  Serial.println(volts, 3);
  delay(10);                               // about 100 samples per second
}
```

## Good data habits

- Record **timestamps** with every sample, not just a sample count.
- Save **metadata**: sensor model, settings, sample rate, units, date, who and where.
- Keep **raw data** untouched; write processed data to new files.
- Use plain, documented formats (CSV for small data, HDF5 for large arrays).

## Real labs

Lab DAQ hardware, high-voltage detector supplies, lasers and cryogenics need institutional training and supervision. At home, stick to phones, simulators and battery-powered kits.

## Common mistakes

- **Sampling too slowly** and analysing an aliased signal.
- **Trusting raw counts** without converting units and calibrating.
- **Losing timing** by assuming a perfectly steady sample rate; phones in particular vary.
- **Averaging away a real signal** with too long a window.

## Recap

- Sensor, conditioning, ADC, timing and storage make up the chain.
- Step size = range / 2^n; resolution is not accuracy.
- Sample faster than twice the highest frequency, and filter first.
- Save timestamps, metadata and raw data.

## Try it

Open the **Do** tab to record accelerometer data with phyphox (or a simulated sensor in Wokwi) and analyse its noise and sample rate.


Explain sensors, analog-to-digital conversion and sampling, and record and analyse real data from a phone or a simulated microcontroller.

Instrumentation and DAQ

Test and Measurement Basics - DAQ | Episode 2 - Define Data Acquisition (DAQ)

Phyphox - Physical Phone Experiments

## Why it matters

Modern physics runs on software. Lattice QCD, climate models, galaxy simulations, gravitational-wave searches and detector reconstruction all need more computing than a laptop, and their results are only as trustworthy as the code. Knowing how **high-performance computing (HPC)** works and how to write **reproducible research software** is now a core physicist skill, and research software engineering is a career in its own right.

## What a cluster is

An HPC cluster is many computers (**nodes**) connected by a fast network.

- **Login nodes:** where you connect with `ssh`, edit files and submit jobs. Never run heavy work here.
- **Compute nodes:** each has many CPU cores (often 32 to 128), lots of memory, and sometimes GPUs.
- **Interconnect:** a low-latency network (such as InfiniBand) so nodes can exchange data quickly.
- **Shared parallel file system:** large, fast storage visible from every node, usually with a scratch area that is purged regularly and is not backed up.

## Schedulers and jobs

You don't run programs directly on compute nodes; you ask a **scheduler** (Slurm is the most common) for resources, and it queues your job until they are free.

```bash
#!/bin/bash
#SBATCH --job-name=heat-scan
#SBATCH --nodes=1
#SBATCH --ntasks=32
#SBATCH --time=02:00:00
#SBATCH --mem=64G

module load python/3.12
srun python run_scan.py --config scan.yaml
```

You submit with `sbatch job.sh`, check the queue with `squeue`, and cancel with `scancel`. Asking for realistic time and memory gets your job started sooner. Computing time is a shared, allocated resource, often granted through proposals.

## Kinds of parallelism

- **Embarrassingly parallel:** many independent runs (parameter scans, Monte Carlo seeds). The easiest and most common; use job arrays.
- **Shared memory (OpenMP, threads):** many cores on one node working on the same arrays.
- **Distributed memory (MPI):** processes on many nodes, each owning part of the problem (for example, part of a grid) and exchanging boundary data with messages.
- **GPUs (CUDA, HIP, or Python tools such as CuPy and JAX):** thousands of simple cores, ideal for dense array maths.

## Amdahl's law

If a fraction `p` of the work can run in parallel on `N` cores, the speed-up is limited:

```text
S = 1 / ((1 - p) + p / N)
p = 0.9, N = 8:    S = 1 / (0.1 + 0.1125) = 4.7
p = 0.9, N = 1000: S = 1 / (0.1 + 0.0009) = 9.9   (never more than 10)
```

The serial part dominates quickly. **Strong scaling** asks how fast a fixed problem gets with more cores; **weak scaling** asks whether a problem that grows with the core count keeps the same run time. Measure both before asking for a large allocation.

## Reproducible research software

A result is reproducible if someone else (or you, in a year) can regenerate it from the same data and code.

- **Version control:** every analysis lives in git; record the commit hash with every result.
- **Pinned environments:** a `requirements.txt` with versions, a conda `environment.yml`, or a container (Apptainer, Docker).
- **Seeds and configuration:** random seeds and parameters live in config files, not in your head.
- **Tests:** check the code against known answers (exact solutions, conservation laws), like the checks in this course.
- **Data provenance:** keep raw data read-only and record each processing step.
- **Sharing and credit:** publish code with a licence, archive a release with a DOI (Zenodo is free), add a `CITATION.cff` file, and follow the FAIR principles (findable, accessible, interoperable, reusable).

## Careers and requirements

National labs, universities and companies employ physicists and research software engineers (RSEs) for this work. Physicist roles need at least a bachelor's degree in physics, and most research positions require a PhD. Many universities offer free HPC training; the Carpentries (https://carpentries.org) teach the core software habits. Checked on 2026-10-23.

## Recap

- Clusters have login and compute nodes, a fast network and shared storage; a scheduler such as Slurm runs your jobs.
- Parallelism ranges from independent runs to MPI and GPUs; Amdahl's law caps the speed-up.
- Version control, pinned environments, seeds, tests and DOIs make research code reproducible.


Explain how clusters, schedulers and parallel code work, and the habits that make research software reproducible.

HPC and research software

High Performance Computing (HPC) - Computerphile

Research Software Engineers and their role for open and reproducible research

Fitting with uncertainties and chi-squared, recording sensor data, and the clusters and habits behind reproducible research code.

Computational & Modern Physics

About this course

Course content

Modern Physics

Computational Methods

Physics Data

Where this course leads

Careers that use it