# Simulation Validation > Validate numerical simulations against closed-form results so a visualisation can be trusted - choosing integrators, comparing vector fields correctly, testing limits by convergence rather than point values, finite-difference step scans, and the stiff-ODE stability limit that turns RK4 into 1e245. - Skill: `bizpers11991-code/simulation-validation` (Agent Skill) - Install (CLI): `npx skillmds@latest add bizpers11991-code/simulation-validation` - Raw SKILL.md: https://api.skillmd.com/api/skills/bizpers11991-code/simulation-validation/raw - Safety review: pending - Works with: Claude Code, Claude.ai, OpenAI Codex - Category: Coding & Dev Tools - Author: bizpers11991-code (https://skillmd.com/u/bizpers11991-code) - Updated: 2026-09-22 - Page: https://skillmd.com/skills/bizpers11991-code/simulation-validation --- # Validating a simulation you intend to show people If you are visualising physics for an audience, **correctness is the product**. This is the pattern that makes that claim checkable instead of aspirational: every solver ships checks against an exact result, all of them report into one registry, and nothing renders until they pass. ## Structure ``` physics/harness.py CHECKS list + check_scalar / check_vectors / record physics/checks_.py one module per solver, importing the harness physics/validate.py imports every check module, prints, exits non-zero ``` The harness lives in its own module so check modules and the runner can both import it without a cycle. Each check appends `(name, ok, note)` — the *note* is what makes a passing suite useful, because it tells you the margin, not just the verdict. ## Compare vector fields by norm, never component-wise ```python err = max(norm(got - want, axis=-1) / maximum(norm(want, axis=-1), 1e-30)) ``` A component that is analytically zero makes a component-wise relative error meaningless: it divides noise by noise and reports a relative error of 1e8 on a field that is perfectly correct. This wastes a lot of time if you meet it first as a mysterious failure. ## Test limits by convergence, not by a point value Some quantities are *first-order results*. They are not supposed to match at any finite parameter value, and a test claiming otherwise is lying about what was verified. Assert that the ratio approaches 1 as the small parameter shrinks, and put the sequence in the note: ``` PASS perihelion advance -> 6piGM/(c^2 a(1-e^2)) measured/Einstein 1.0081 -> 1.0020 -> 1.0005 as c doubles PASS capacitor E -> sigma/eps0 as the gap closes 0.8623 -> 0.9229 -> 0.9555 at gap/size 0.20, 0.10, 0.05 PASS ring sum -> closed-form drag k_sum/k_closed 0.9386 -> 1.0000 -> 1.0000 ``` Each of those started as a failing point-comparison. Each time, the simulation was right and the *test* was wrong about what it should equal. Before loosening a tolerance, work out whether the quantity is a limit. Errors falling by 4x per halving of the small parameter is the signature of a first-order expansion — good confirmation you understand the discrepancy rather than papering over it. ## Choose the integrator for the invariant you care about | Problem | Use | Why | |---|---|---| | charged particle in B | **Boris pusher** | conserves energy; RK4/Euler spiral outward and the visual is simply wrong | | gravitational orbits | **velocity Verlet** | symplectic; energy oscillates in a bounded band instead of drifting | | velocity-dependent forces (Coriolis) | RK4 | Verlet cannot absorb them; police it with an exact invariant instead | | field-line tracing | RK4 on the **normalised** field | parameter becomes arc length, so step size is a real distance | Measured: Boris held `|v|` to 1.8e-15 over 1500 steps; Verlet's orbital energy drifted 1.5e-14 over 60 orbits. RK4 on the same orbit loses energy monotonically and spirals in — at 30 fps that looks like physics and is numerics. When you must use RK4, check a conserved quantity it does not know about: the Jacobi constant for the restricted three-body problem held to 3e-14. ## Stiff ODEs: RK4 explodes, it does not wobble RK4 on `dv/dt = g - v/tau` is unstable beyond about `dt = 2.78 tau`. Past that it does not degrade gracefully — it returns **1e245** and keeps going. A renderer will happily build a scene from that and produce an empty frame. This is easy to hit when the settle time is milliseconds and "one frame" is tens of them. Make the integrator refuse: ```python tau = mass / k if dt > 0.5 * tau: raise ValueError(f"dt={dt:.4g} exceeds half the settle time tau={tau:.4g}; " "RK4 is unstable here. Integrate finely and resample.") ``` Integrate at `tau/20` and `np.interp` onto your keyframe times. ## Finite differences need a step scan Checking something like `div B = 0` numerically has no single correct `h`: too large and truncation error dominates, too small and cancellation does. Scan `h` over several decades and take the best, reporting it as such: ``` PASS div B = 0 (dipole) best |div B| / max|dBi/dxi| = 2.50e-08 over h scan ``` Normalise by the largest individual derivative, or you are testing the field's magnitude rather than its divergence. ## Fit, do not average Recovering an orbit radius by taking the mean of a *partial* orbit gives the wrong centre. Use a least-squares circle fit, which is robust to partial arcs: ```python A = column_stack([2*x, 2*y, ones(len(x))]) cx, cy, c = lstsq(A, x**2 + y**2)[0] radius = sqrt(c + cx**2 + cy**2) ``` Also sample over whole periods where you can, and compare a chord to the field at the chord's **midpoint** rather than its start — comparing at the start introduces an O(h) bias that looks like a real error. ## Exclude singularities from tracing Every field line of a point dipole converges on its singularity, where a fixed-step integrator produces garbage — and the point-dipole model is invalid inside the magnet's body anyway. Pass the source bodies as `(centre, radius)` exclusion blobs and stop a line when it enters one. ## Check the vectorised path against the naive one Keep a slow, obviously-correct reference implementation (a per-source Python loop) and assert the fast einsum path matches it exactly on random points: ``` PASS vectorised field == reference sum rel err 0.00e+00 over 300 random points ``` This is the cheapest check in the suite and the one most likely to catch a silent indexing mistake during optimisation. ## When a check fails, suspect the check In this build, failing checks broke down roughly as: - **the test was wrong about the expected value** (a limit, or an unphysical parameter choice) — most common - **a genuine sign error** — a `theta-hat` built as `z - cos(th) r` instead of `cos(th) r - z`, which is `-theta-hat`. Watch for a relative error of **exactly 2.0**: that is `got = -want`, and it is a sign flip every time - **a real numerical problem** — the stiff-ODE blow-up above A relative error of exactly 2.0, or a ratio of exactly -1, should send you looking for a sign before anything else. ## Parameters are part of the physics A solver can be perfect and the *scenario* nonsense. Two from this build: - A magnetic bottle with a gyroradius **51% of the coil radius** is not adiabatic; the magnetic moment drifted 380% and the particle was not really trapped. At 6% it conserved to 3.3%. - A copper-tube eddy brake with an implausible bore and a weak magnet gave a terminal velocity of 4.9 m/s. Real parameters — a half-inch pipe and a stack of three 10 mm magnets — gave 5.9 cm/s, matching the lecture demo. Sanity-check derived quantities against something you can look up, and put the number in the check's note so it stays visible. ## Derive rather than quote, then close the loop Where a closed form is needed, derive it in the module docstring and then verify it against an independent numerical route that shares none of the algebra. The eddy-drag coefficient here was derived by integrating the ring kernel, then confirmed by a discrete ring sum converging on it to 1e-4 — and the ring kernel itself was checked against the already-validated Biot-Savart core. A quoted constant you half-remember is a liability. A derivation with a convergence test behind it is not.