Numerics · Viv posts

Skip to content

This site is version 0.0.1. Most of the code on this site is drafted by an agent and needs review. I apologize for the AI-speak, we are working on rewriting it!

Floating point

demoAINumericswritten with an agent#5589c09a31

Adam Zettel@adamSep 25, 2026

Three floating-point failures in summation, cancellation and RANDU

fl (a + b) = (a + b) (1 + δ), ∣ δ ∣ \leq 2^{- 53}; x_{k + 2} = 6 x_{k + 1} - 9 x_{k} mod 2^{31}

Adding $1/ k$ for $k$ up to a million left to right ends 414 units in the last place (ulps) from the correctly rounded sum; Kahan's compensated sum is correctly rounded at all 31 checkpoints from 1 to $1 0^{6}$ . Evaluating $(x - 1)^{7}$ from its expanded coefficients near $x = 1$ gives rounding noise: 89 sign changes among 401 points, where the polynomial has one, and 8.5 times the true value at $x = 1.005$ . The factored form is within $1.4 ε$ everywhere. The random number generator RANDU puts every triple of successive draws on one of 15 planes.

Sum of 1/k, k = 1 .. N: error in ulps of the exact sum

left to rightKahan

$(x - 1)^7$: factored, and expanded by Horner's rule

0.98 to 1.02

Horner, expanded

0.995 to 1.005

Horner, expanded

$9u_k - 6u_{k+1} + u_{k+2}$ against $u_k$, 600 triples

RANDU

MINSTD

✓ 15 tests passchecked when the site was built#floating-point #rounding #kahan-summation #horner #cancellation #random-numbers #notebook

ReadOpen in the workbench

Look inside this program

Atlas
The program’s map: every function and which calls which.
There’s no recorded run of this program yet.

Integrating in time

demoAINumericswritten with an agent#70ee9cc956

Adam Zettel@adamSep 21, 2026

Forward Euler loses the planet

r ¨ = - \frac{r}{∣ r ∣ ^{3}}, E = \frac{1}{2} ∣ r ˙ ∣^{2} - \frac{1}{∣ r ∣}

Ten orbits of a planet on an ellipse, 200 steps an orbit, stepped two ways at the same cost. Forward Euler gains energy until the planet is unbound. Leapfrog's energy error stays under one percent, and it keeps the angular momentum to rounding.

Energy error over ten orbits, equal step size

forward Eulerleapfrog

✓ 7 tests passchecked when the site was built#orbits #integrators

ReadOpen in the workbench

Look inside this program

Atlas
The program’s map: every function and which calls which.
There’s no recorded run of this program yet.

demoAINumericswritten with an agent#82acda210a

Adam Zettel@adamSep 25, 2026

Forward and backward Euler on either side of h = 2/λ

y^{'} = - λ (y - cos t) - sin t, ∣1 - hλ ∣ < 1 ⟺ h < 2/ λ

The exact solution of this equation is $cos t$ , which a step of 0.1 would follow to a few digits. With $λ = 1000$ , forward Euler is within $4 \times 1 0^{- 7}$ of it at $t = 2$ with $h = 0.0019$ , and off by $2.7 \times 1 0^{33}$ at $h = 0.0021$ . A disturbance is multiplied by $1 - hλ$ each step, and the measured factor matches $∣1 - hλ ∣$ to $2 \times 1 0^{- 14}$ at 36 step sizes. Backward Euler's factor is $1/ (1 + hλ)$ , below 1 at every step size, and at $h = 0.1$ it is still within $2 \times 1 0^{- 5}$ of $cos 2$ .

cos t, exactforward Euler, h = 0.0019forward Euler, h = 0.0021backward Euler, h = 0.01

forward Euler, measuredbackward Euler, measured

✓ 9 tests passchecked when the site was built#ode #stiffness #forward-euler #backward-euler #stability #notebook

ReadOpen in the workbench

Look inside this program

Atlas
The program’s map: every function and which calls which.
There’s no recorded run of this program yet.

problemAINumericswritten with an agent#e6f3a56fdf

Adam Zettel@adamSep 25, 2026

Make Euler stable on a stiff equation

y^{'} = - 50 y, h = 0.05 : y_{n + 1} = (1 - 50 h) y_{n} = - 1.5 y_{n}

The solution of $y^{'} = - 50 y$ decays to zero, but forward Euler with step $h = 0.05$ multiplies $y$ by $1 - 50 h = - 1.5$ every step, and after 12 steps it is 129.7 instead of $1 0^{- 13}$ . Rewrite decay_step and system_step as backward Euler, which takes the slope at the end of the step instead of the start, so that no step size makes them grow.

one step of 0.05 multiplies y by 1 / (1 + 2.5) = 0.285714
twelve steps of 0.05 leave 2.9593e-7: small and positive, never negative
a step 25 times the forward-Euler limit, h = 1, gives 1/51 = 0.0196
ten steps of 0.1 on the system end at numpy's (0.38554, 0.38554): the fast mode is gone

✗ 4 tests to solvethe unfinished functions run, and the failing tests show what they must computeno solutions yet #ode #stiffness #implicit-methods #backward-euler #problem

ReadOpen in the workbench

Look inside this program

Atlas
The program’s map: every function and which calls which.
There’s no recorded run of this program yet.

problemAINumericswritten with an agent#dfc71ed5d9

Adam Zettel@adamSep 25, 2026

Make Numerov find the decaying solution

y^{''} = y, y (x) = e^{- x} : outward, the error grows as e^{2 x}

Numerov's method integrates $y^{''} = y$ outward from $y (0) = 1$ , $y (h) = e^{- h}$ , aiming at the decaying solution $e^{- x}$ . By $x = 9$ it is off by a factor of 30, and by $x = 18$ by $2 \times 1 0^{9}$ . Rewrite decaying so that it follows $e^{- x}$ to Numerov's own accuracy, within $2.7 \times 1 0^{- 5}$ , on both intervals.

with h = 0.15 on [0, 9] every point is within 1.3e-5 of e^(-x), where numpy's inward run reaches 8.87e-6
with h = 0.075 on [0, 9] every point is within 8e-7, where numpy reaches 5.58e-7
with h = 0.15 on [0, 18], where outward is off by 2e9, every point is within 2.7e-5 (numpy reaches 1.84e-5)
y(0) = 1
halving h divides the error by 12 to 20, as a fourth-order method's should (numpy: 15.9)

✗ 3 tests to solve2 pass already · the unfinished functions run, and the failing tests show what they must computeno solutions yet #ode #numerov #shooting #stability #problem

ReadOpen in the workbench

Look inside this program

Atlas
The program’s map: every function and which calls which.
There’s no recorded run of this program yet.

Approximation

articleAINumericswritten with an agent#b93d4ae27a

Adam Zettel@adamSep 21, 20261 min read

The Runge phenomenon, where more points give a worse fit

f (x) = \frac{1}{1 + 25 x ^{2}}, ∥ f - p_{n} ∥_{\infty} \to \infty

Interpolate $1/ (1 + 25 x^{2})$ by a polynomial through equally spaced points and the error grows as you add points. Move the same number of points to $cos (j π / n)$ and it falls geometrically, at a rate you can work out from where the function's poles are.

equally spacedChebyshev

1 paragraph · 4 functions · 1 min

✓ 3 tests passchecked when the site was built#interpolation #approximation

Read the articleOpen in the workbench

Look inside this program

Atlas
The program’s map: every function and which calls which.
There’s no recorded run of this program yet.

articleAINumericswritten with an agent#a67916e929

Adam Zettel@adamSep 21, 20262 min read

The Gibbs overshoot never shrinks

x max ∣ S_{m} (x) - f (x) ∣ \to \frac{2}{π} \int_{0}^{π} \frac{sin t}{t} d t - 1

Add more terms to the Fourier series of a square wave and the wiggle at the jump gets narrower, never shorter. It overshoots by 8.95% of the jump however many terms you take, and the number comes from an integral.

Square wave at the jump, partial sums m = 5 to 80

0.0 / 12.5 s

8 paragraphs · 3 functions · 2 min

✓ 4 tests passchecked when the site was built#fourier #analysis

Read the articleOpen in the workbench

Look inside this program

Atlas
The program’s map: every function and which calls which.
There’s no recorded run of this program yet.

Quadrature

demoAINumericswritten with an agent#0b475cd75c

Adam Zettel@adamSep 25, 2026

Gauss-Legendre quadrature against tanh-sinh on a smooth and a singular integral

\int_{- 1}^{1} p (x) d x = i = 1 \sum n w_{i} p (x_{i}) for every polynomial p of degree \leq 2 n - 1

Five Gauss-Legendre points integrate every polynomial of degree 9 or less exactly, and miss degree 10 by exactly what the error formula says. On $e^{x}$ seven points reach rounding. On $1/ 1 - x^{2}$ , which is infinite at both ends of the interval, 128 points still miss by $1.4 \times 1 0^{- 2}$ , and the error only halves when the points double. The tanh-sinh rule gets that integral to rounding with 37 points.

5 points3 points

Gauss-Legendretanh-sinh

Gauss-Legendretanh-sinh

✓ 13 tests passchecked when the site was built#quadrature #integration #gauss-legendre #tanh-sinh #notebook

ReadOpen in the workbench

Look inside this program

Atlas
The program’s map: every function and which calls which.
There’s no recorded run of this program yet.

Eigenvalues

demoAINumericswritten with an agent#8b2f1cac21

Adam Zettel@adamSep 21, 2026

How far a finite-difference Hamiltonian's levels sit below n + 1/2

H^= - \frac{1}{2} \partial_{x}^{2} + \frac{1}{2} x^{2}, E_{n} = n + \frac{1}{2}

Put $H = - \frac{1}{2} \partial_{x}^{2} + \frac{1}{2} x^{2}$ on a grid and diagonalize: the levels come out at $n + \frac{1}{2}$ , a little low. How low is predictable. First-order perturbation theory gives $- \frac{h ^{2}}{32} (2 n^{2} + 2 n + 1)$ , and for the six lowest levels the measured shift matches it to within half a percent (the test allows 4%).

The three lowest eigenfunctions

✓ 3 tests passchecked when the site was built#quantum #eigenvalues

ReadOpen in the workbench

Look inside this program

Atlas
The program’s map: every function and which calls which.
There’s no recorded run of this program yet.

demoAINumericswritten with an agent#fd96a8df95

Adam Zettel@adamSep 21, 2026

The eigenvalues of a random matrix fill Wigner's semicircle

ρ (λ) = \frac{1}{2 π σ ^{2}} 4 σ^{2} - λ^{2}, ∣ λ ∣ \leq 2 σ

Fill a symmetric $200 \times 200$ matrix with independent random numbers, divide by $N$ , and look at its eigenvalues. They do not scatter. They fill the semicircle $\frac{1}{2 π} 4 - x^{2}$ , whatever distribution the entries came from.

Eigenvalues of one N x N random symmetric matrix

histogramsemicircle

✓ 4 tests passchecked when the site was built#random-matrices #statistics

ReadOpen in the workbench

Look inside this program

Atlas
The program’s map: every function and which calls which.
There’s no recorded run of this program yet.

articleAINumericswritten with an agent#e20fe3d0a4

Adam Zettel@adamSep 24, 20265 min read

LOBPCG finds the lowest states of a particle in a box

λ_{1} = x \neq = 0 min \frac{x ^{H} H x}{x ^{H} x}

The six lowest states of a particle in a square box, from random starting vectors, in 103 steps of LOBPCG, and the 576 by 576 Hamiltonian is never written down. Each step searches three blocks of vectors: the current guesses, their residuals, and the step each guess took last time. Leave out the third block and the same run takes 797 steps.

Residual norm of each state, per LOBPCG step

(1,1)(1,2)(2,1)(2,2)(1,3)(3,1)tolerance

The worst wanted residual, with and without the previous step p

LOBPCG: x, w and p (held at its last value once done)steepest descent: x and w

The six lowest states

(1,1)

24 × 24 grid

the (1,2) and (2,1) pair

24 × 24 grid

the (1,2) and (2,1) pair

24 × 24 grid

(2,2)

24 × 24 grid

the (1,3) and (3,1) pair

24 × 24 grid

the (1,3) and (3,1) pair

24 × 24 grid

26 paragraphs · 28 functions · 5 min

✓ 6 tests passchecked when the site was built#eigenvalues #linear-algebra #quantum

Read the articleOpen in the workbench

Look inside this program

Atlas
The program’s map: every function and which calls which.
There’s no recorded run of this program yet.

demoAINumericswritten with an agent#295447f910

Adam Zettel@adamSep 25, 2026

Jacobi's eigenvalue method

off (A) = (i \neq = j \sum a_{ij}^{2})^{1/2}

Jacobi's method diagonalizes a symmetric matrix by plane rotations, each of which zeroes one off-diagonal pair. On an 8 by 8 matrix the off-diagonal norm goes 1.33, 0.453, 0.0299, $3.2 \times 1 0^{- 5}$ , $8 \times 1 0^{- 13}$ over four sweeps, and after the fifth it is under $1 0^{- 15}$ , rounding error. The diagonal is then the spectrum, and it matches numpy.linalg.eigh to $1 0^{- 14}$ .

Off-diagonal norm after each Jacobi sweep, 8 by 8 matrix

off(A)

✓ 10 tests passchecked when the site was built#eigenvalues #linear-algebra #jacobi #convergence #notebook

ReadOpen in the workbench

Look inside this program

Atlas
The program’s map: every function and which calls which.
There’s no recorded run of this program yet.

problemAINumericswritten with an agent#9fbf3ac156

Adam Zettel@adamSep 25, 2026

Make power iteration converge on a tied pair

A = (0220), λ = \pm 2 : v_{k + 1} = \frac{( A + s I ) v _{k}}{∥ ( A + s I ) v _{k} ∥}

Power iteration multiplies a vector by $A$ again and again and finds the eigenvalue of largest magnitude. The matrix above has two, $+ 2$ and $- 2$ , of the same magnitude, and the iterate swaps its two entries every step forever. Write shifted_power so that the shift $s$ breaks the tie: iterate with $A + s I$ , whose eigenvalues are $λ + s$ .

shift +1 finds eigenvalue 2, vector (1, 1) / sqrt 2
shift -1 finds eigenvalue -2, vector (1, -1) / sqrt 2
on [[1, 2], [2, -1]] shift +1 finds sqrt 5 = 2.2360680 with numpy's vector
and shift -1 finds -sqrt 5

✗ 4 tests to solvethe unfinished functions run, and the failing tests show what they must computeno solutions yet #eigenvalues #power-iteration #shifts #linear-algebra #problem

ReadOpen in the workbench

Look inside this program

Atlas
The program’s map: every function and which calls which.
There’s no recorded run of this program yet.

demoAINumericswritten with an agent#d003d94d29

Adam Zettel@adamSep 25, 2026

Lanczos finds the largest eigenvalue three times

A q_{j} = β_{j - 1} q_{j - 1} + α_{j} q_{j} + β_{j} q_{j + 1}

Sixty steps of Lanczos on a 200 by 200 symmetric matrix whose eigenvalues are known exactly. Both ends of the spectrum, $\pm 4$ , come out right to $1 0^{- 14}$ . But in floating point the basis stops being orthogonal once the top eigenvalue has converged, the direction of its eigenvector creeps back in, and the eigenvalue 4, which the matrix has once, is found again and again. Orthogonalizing each new vector against all the old ones stops it.

Plain Lanczos against full reorthogonalization, step by step

Component of the new basis vector along the eigenvector for 4

plain three-term recurrencefull reorthogonalization

Largest overlap of the new basis vector with an earlier one

plain three-term recurrencefull reorthogonalization

The three largest Ritz values

plain, largestplain, 2ndplain, 3rdreorthogonalized, largestreorthogonalized, 2ndreorthogonalized, 3rd

✓ 12 tests passchecked when the site was built#eigenvalues #linear-algebra #lanczos #krylov #rounding #notebook

ReadOpen in the workbench

Look inside this program

Atlas
The program’s map: every function and which calls which.
There’s no recorded run of this program yet.

Orthogonalization

problemAINumericswritten with an agent#386cd64c28

Adam Zettel@adamSep 25, 2026

Make Gram-Schmidt orthogonal to rounding

a_{k} = e_{0} + ε e_{k + 1}, ε = 1 0^{- 8} : i \neq = j max ∣ q_{i} \cdot q_{j} ∣ = 0.5 (classical)

Ten nearly parallel vectors, each $e_{0}$ plus a tiny $1 0^{- 8}$ in its own direction. Classical Gram-Schmidt turns them into ten "orthonormal" vectors of which two are at 60 degrees: their dot product is 0.5. Rewrite orthonormalize so the vectors it returns are orthogonal to $1 0^{- 15}$ and equal the exact ones to $1 0^{- 12}$ .

the rows are orthonormal to 1e-15 (classical Gram-Schmidt gives 0.5, modified 7.1e-9)
the second vector is mpmath's (7.07e-9, -1/sqrt 2, 1/sqrt 2, 0, ...)
the last vector is mpmath's (1.05e-9, -0.10541 nine times, 0.94868)
each input row is rebuilt from the rows before it: a_k = sum over j <= k of (q_j . a_k) q_j

✗ 4 tests to solvethe unfinished functions run, and the failing tests show what they must computeno solutions yet #linear-algebra #gram-schmidt #orthogonalization #rounding #problem

ReadOpen in the workbench

Look inside this program

Atlas
The program’s map: every function and which calls which.
There’s no recorded run of this program yet.

Partial differential equations

problemAINumericswritten with an agent#e1e71349a4

Adam Zettel@adamSep 25, 2026

Make the heat-equation step stable past r = 1/2

u_{i}^{n + 1} = u_{i}^{n} + r (u_{i - 1}^{n} - 2 u_{i}^{n} + u_{i + 1}^{n}), r = \frac{Δ t}{Δ x ^{2}} = 0.6 \Rightarrow blows up

The explicit (FTCS) step for $u_{t} = u_{xx}$ multiplies the grid's most jagged mode by about $1 - 4 r$ each step, and for $r = 0.6$ that is $- 1.4$ : a checkerboard that should fade grows to 215 in 16 steps. Rewrite heat_step as backward Euler, which needs a tridiagonal solve per step and damps every mode for every $r$ .

the smoothest mode decays by 1 / (1 + 2.4 sin^2(pi/40))^16 = 0.790844 at r = 0.6
the most jagged mode, m = 19, decays to 3.362e-9 of itself instead of growing
the checkerboard at r = 0.6 fades to numpy's 0.0068342 after 16 steps
and at r = 5, ten times the FTCS limit, to numpy's 0.0012278
the ends stay at 0

✗ 4 tests to solve1 pass already · the unfinished functions run, and the failing tests show what they must computeno solutions yet #pde #heat-equation #stability #implicit-methods #problem

ReadOpen in the workbench

Look inside this program

Atlas
The program’s map: every function and which calls which.
There’s no recorded run of this program yet.

Matrix equations

demoAINumericswritten with an agent#f45cc94837

Adam Zettel@adamSep 25, 2026

The error of a linear solve follows the condition number

\frac{∥ x ^ - x ∥}{∥ x ∥} ≲ κ (A) ε, κ (A^{T} A) = κ (A)^{2}

Gaussian elimination with partial pivoting solves Hilbert systems whose exact solution is all ones. From $n = 3$ to 12 the condition number $κ$ grows from 524 to $1.7 \times 1 0^{16}$ , and the error grows with it, staying between 0.005 and 0.18 times $κ ε$ , while the residual stays under $ε$ . In least squares, fitting a degree-11 polynomial to 40 points, Householder QR is off by $3.6 \times 1 0^{- 9}$ (about $κ ε$ ), and the normal equations, which square $κ$ , by 0.18.

GEPP errorresidual / (|A| |x|)

Householder QRnormal equations

✓ 11 tests passchecked when the site was built#linear-algebra #conditioning #hilbert-matrix #least-squares #householder-qr #normal-equations #notebook

ReadOpen in the workbench

Look inside this program

Atlas
The program’s map: every function and which calls which.
There’s no recorded run of this program yet.

demoAINumericswritten with an agent#cec2f50a64

Adam Zettel@adamSep 25, 2026

Jacobi, SOR and conjugate gradients on the Poisson equation

ρ_{Jacobi} = cos π h, ρ_{SOR} = \frac{1 - sin π h}{1 + sin π h}, κ = cot^{2} \frac{π h}{2}

The five-point Poisson equation on a 15 by 15 grid, solved to a relative residual of $1 0^{- 8}$ four ways. Jacobi takes 948 sweeps, Gauss-Seidel 475, SOR with the optimal factor 60, and conjugate gradients 25 iterations. Each count follows from one number the theory predicts: Jacobi's error shrinks by $cos π h$ per sweep once the smoothest mode dominates, Gauss-Seidel's by its square, and SOR and CG depend on $κ \propto 1/ h$ . Halve $h$ and Jacobi takes four times as many sweeps; SOR and CG take twice as many.

JacobiGauss-SeidelSORCG

Iterations to a residual of 1e-8 against grid size

JacobiSORCG

✓ 15 tests passchecked when the site was built#linear-algebra #iterative-methods #poisson #sor #conjugate-gradient #notebook

ReadOpen in the workbench

Look inside this program

Atlas
The program’s map: every function and which calls which.
There’s no recorded run of this program yet.

demoAINumericswritten with an agent#7e2c4a61ad

Adam Zettel@adamSep 23, 2026

Solving the Riccati equation without integrating it

A^{⊤} P + P A - P B R^{- 1} B^{⊤} P + W = 0

Every LQR controller comes from a quadratic matrix equation, and the tempting way to solve it is to integrate the differential version backwards until it stops moving. On a one-state system that method returns 73.783 where the answer is 414.214, with no warning. The Schur method Laub published in 1979 returns the answer to rounding error, and this post runs it.

the six stable ones, which P is built fromthe six unstable ones, discarded

Each subdiagonal entry shrinking to zero

entry 9, the first to goentry 7, which goes nextentry 3, higher up

The subdiagonal being driven to zero

0.0 / 9.9 s

Relative residual, ten systems, two methods

the Schur methodthe matrix sign function

✓ 10 tests passchecked when the site was built#control #linear-algebra #schur #riccati

ReadOpen in the workbench

Look inside this program

Atlas
The program’s map: every function and which calls which.
There’s no recorded run of this program yet.

PDEs

demoAINumericswritten with an agent#a4811d5668

Adam Zettel@adamSep 25, 2026

Upwind advection smears the pulse by its modified equation

u_{t} + c u_{x} = D u_{xx}, D = \frac{c h}{2} (1 - ν)

First-order upwind carries a Gaussian pulse once around a periodic domain at Courant number $ν = 0.5$ . The mass is kept to rounding and no new maximum or minimum appears, but the peak falls from 0.9975 to 0.5767. The scheme's modified equation adds a diffusion term $D u_{xx}$ with $D = c h (1 - ν) /2$ , which predicts the fall to $1/ 3 = 0.5774$ and predicts the pulse's variance growth, $2 D T$ , exactly.

start, and the exact answer after one transitupwind after 400 stepsthe modified equation's Gaussian

Upwind after one transit, against its modified equation

Peak height

upwind

Variance growth

upwind, measured

✓ 10 tests passchecked when the site was built#pde #advection #upwind #numerical-diffusion #modified-equation #notebook

ReadOpen in the workbench

Look inside this program

Atlas
The program’s map: every function and which calls which.
There’s no recorded run of this program yet.

demoAINumericswritten with an agent#9458e4aca8

Adam Zettel@adamSep 25, 2026

Crank-Nicolson is second order in time, backward Euler first

\frac{u ^{n + 1} - u ^{n}}{Δ t} = θ δ^{2} u^{n + 1} + (1 - θ) δ^{2} u^{n}, θ = \frac{1}{2} or 1

The heat equation $u_{t} = u_{xx}$ from $u (x, 0) = x (1 - x)$ , stepped to $t = 0.1$ with six time steps from 0.02 down to 0.000625. Each halving cuts the Crank-Nicolson error by a factor of 4.0 to 4.7 (measured order 2.09 to 2.23, then 2.000) and the backward Euler error by 1.9 to 2.0 (order 0.94 to 0.996). The errors are measured against the exact solution of the space-discretized equation, a discrete sine series. Against the exact Fourier series of the PDE, the Crank-Nicolson error stops falling at the space error, $1.95 \times 1 0^{- 6}$ .

Error at t = 0.1 against the time step

backward Eulerslope 1Crank-Nicolsonslope 2Crank-Nicolson against the PDE series

✓ 9 tests passchecked when the site was built#pde #heat-equation #crank-nicolson #backward-euler #order-of-accuracy #notebook

ReadOpen in the workbench

Look inside this program

Atlas
The program’s map: every function and which calls which.
There’s no recorded run of this program yet.

demoAINumericswritten with an agent#7389123e4c

Adam Zettel@adamSep 25, 2026

Burgers' equation in nonconservative form puts the shock in the wrong place

s = \frac{f ( u _{L} ) - f ( u _{R} )}{u _{L} - u _{R}} = \frac{u _{L} + u _{R}}{2}, f (u) = \frac{1}{2} u^{2}

For smooth solutions, $u_{t} + (u^{2} /2)_{x} = 0$ and $u_{t} + u u_{x} = 0$ are the same equation. For a shock from $u_{L} = 1$ to $u_{R} = 0.25$ they are not. A scheme built on the first form, Godunov's, moves the shock at 0.625, the Rankine-Hugoniot speed $(u_{L} + u_{R}) /2$ , to rounding. Upwinding the second form moves it at 0.5132, on 400, 800 and 1600 cells alike, so refining the grid does not fix it. With $u_{R} = 0$ the nonconservative shock stays at $x = 0$ , with speed 0.

start, t = 0exact: shock at x = 0.625conservative (Godunov)

Rankine-Hugoniot, x = 0.625 tconservativenonconservative

✓ 7 tests passchecked when the site was built#pde #burgers #shocks #conservation-laws #rankine-hugoniot #godunov #notebook

ReadOpen in the workbench

Look inside this program

Atlas
The program’s map: every function and which calls which.
There’s no recorded run of this program yet.

Roots

demoAINumericswritten with an agent#cd6c853aa2

Adam Zettel@adamSep 21, 2026

Which root does Newton's method find?

z_{k + 1} = z_{k} - \frac{z _{k}^{3} - 1}{3 z _{k}^{2}}

Solve $z^{3} = 1$ by Newton's method from every point of a square and color each start by the root it reaches. Near a root the answer is obvious. Along the lines between them it is a fractal: every boundary point touches all three basins.

Basins of the three cube roots of 1

Im z

Re z

128 × 128 grid

✓ 3 tests passchecked when the site was built#fractals #root-finding

ReadOpen in the workbench

Look inside this program

Atlas
The program’s map: every function and which calls which.
There’s no recorded run of this program yet.

Also in Numerics

articleAINumericswritten with an agent#1173edd4ea

Adam Zettel@adamOct 14, 20262 min read

Finite differences and spectral methods

f^{'} (x_{j}) \approx \frac{f _{j + 1} - f _{j - 1}}{2 h} vs f^{'} (x_{j}) \approx k \sum ik f^_{k} e^{ik x_{j}}

Two ways to differentiate a function you only know at grid points. A central difference uses the two neighbors; a spectral derivative uses every point. On a smooth periodic function the spectral one reaches rounding error with a few dozen points. Put a jump in the function and, away from the jump, the central difference is the accurate one.

Largest derivative error against the number of points

spectral, exp(sin x)central difference, exp(sin x)spectral, sin x + step, far pointscentral difference, sin x + step, far points

10 paragraphs · 8 functions · 2 min

✓ 3 tests passchecked when the site was built#numerics #fourier #derivatives

Read the articleOpen in the workbench

Look inside this program

Atlas
The program’s map: every function and which calls which.
There’s no recorded run of this program yet.

demoAINumericswritten with an agent#ca6b837769

Adam Zettel@adamSep 25, 2026

Electron mobility fitted to NIST's certified digits

b min i \sum (\frac{b _{1} + b _{2} x _{i} + b _{3} x _{i}^{2} + b _{4} x _{i}^{3}}{1 + b _{5} x _{i} + b _{6} x _{i}^{2} + b _{7} x _{i}^{3}} - y_{i})^{2}

A seven-parameter rational model fitted to NIST's semiconductor electron-mobility data by lmder, the Levenberg-Marquardt routine of MINPACK, translated to Viv step by step. Every trial step matches a gfortran build of the original Fortran to the last bit, the answer matches scipy.optimize.leastsq, and the parameters land within 0.001 standard deviations of NIST's certified values.

NIST's measurementsthe fitted rational model

Convergence, one point per trial step

How far the norm of the residuals is from its final value

norm of residuals - final norm

The trust radius

✓ 8 tests passchecked when the site was built#least-squares #levenberg-marquardt #minpack #curve-fitting #nist

ReadOpen in the workbench

Look inside this program

Atlas
The program’s map: every function and which calls which.
Explore
Run it and step through what each part did, with the code, values and plots beside it.

demoAINumericswritten with an agent#009c306ada

Adam Zettel@adamSep 24, 2026

Newton's method doubles the correct digits

e_{n + 1} \approx \frac{e _{n}^{2}}{2 2}

Newton's method on $x^{2} = 2$ , one cell at a time. From $x_{0} = 1$ the error is $1.6 \times 1 0^{- 12}$ after four steps and zero after five; bisection needs 39 halvings to promise as much. Each Newton step squares the error, and a cell checks the constant in front.

bisection, half the interval

✓ 7 tests passchecked when the site was built#newton #root-finding #convergence #notebook

ReadOpen in the workbench

Look inside this program

Atlas
The program’s map: every function and which calls which.
There’s no recorded run of this program yet.