PDE solvers
Partial differential equations, solved and checked
15 equations. Each is one Viv program that solves it, draws the solution, and ends in tests against something exact: a closed-form solution, a conserved quantity, a measured order of accuracy. 57 tests in all.
The programs are short enough to read whole, and every page has the source in an editable cell: change the scheme, the time step or the initial condition and run it again.
Diffusion (parabolic)
Smoothing, irreversible, and stiff: an explicit step is limited by , so the working methods are implicit and the work is a tridiagonal solve.
The heat equation
Crank-Nicolson: implicit, second order in space and time, stable at any time step
checked: the exact decay of a sine mode; the measured order of accuracy; and that the explicit scheme blows up exactly where theory says it must
3 tests · 82 lines
The heat equation in two dimensions
alternating-direction implicit (Peaceman-Rachford): one tridiagonal solve per grid line, per half-step
checked: two decaying modes against their exact rates, the measured order, and a time step 128 times the explicit limit
3 tests · 101 lines
The Black-Scholes equation
Crank-Nicolson on a grid in the stock price, started with four implicit half-steps (Rannacher)
checked: a European call against the closed form: the price, the delta, and the measured order with and without the start-up steps
4 tests · 126 lines
Waves and transport (hyperbolic)
Information moves at a finite speed along characteristics. Nothing smooths the solution, so the scheme's own dissipation and dispersion are what you see, and shocks need conservation.
The advection equation
upwind (first order) and Lax-Wendroff (second order), on a periodic grid
checked: that one trip around the ring returns the pulse to where it started; and both measured orders
3 tests · 73 lines
The wave equation
leapfrog: centered in time and space, second order, explicit
checked: a standing wave against its exact solution; the measured order; and that the energy does not drift
3 tests · 89 lines
A shock and a rarefaction in inviscid Burgers' equation
Godunov's method: the exact Riemann solution at every cell face
checked: the shock speed, exact conservation, and the rarefaction fan u = x/t
4 tests · 85 lines
The Euler equations of gas dynamics
finite volumes with the HLL approximate Riemann flux, first order; the exact solution is computed alongside
checked: Sod's shock tube: the exact solution against the jump conditions, then the computed plateau, shock position, conservation and L1 convergence against it
7 tests · 245 lines
Steady states and modes (elliptic)
No time to march through: every unknown depends on every other, and the work is a large sparse linear system or its eigenvalues.
Laplace's equation
successive over-relaxation, with the optimal relaxation factor for a square
checked: a value symmetry fixes exactly; and that over-relaxation beats plain Gauss-Seidel by the factor theory predicts
3 tests · 63 lines
Poisson's equation
five-point Laplacian, solved by conjugate gradients without ever forming the matrix
checked: a manufactured solution; the measured order; and that conjugate gradients converges in far fewer steps than unknowns
3 tests · 95 lines
The modes of a drum, a Helmholtz eigenproblem
the five-point Laplacian as a dense symmetric matrix, diagonalized by eig_sym
checked: the discrete eigenvalues against their closed form (to rounding), the continuum ones to second order, and the degeneracy a square must have
3 tests · 82 lines
Waves whose speed depends on wavelength (dispersive)
The linear part is exact in Fourier space for any time step, so these are split-step spectral codes built on fft.
The time-dependent Schrödinger equation
split-step Fourier (Strang splitting): the kinetic step is exact in momentum space
checked: a free packet's spreading, a coherent state's oscillation, norm conservation, and the measured order of the splitting
6 tests · 108 lines
The Korteweg-de Vries equation
Fourier pseudo-spectral in space; Strang splitting in time, with the dispersive step exact
checked: a soliton against its exact solution: shape, speed, conserved mass, and the measured order
5 tests · 163 lines
The nonlinear Schrödinger equation
split-step Fourier: the nonlinear step is ALSO exact, because it conserves the modulus
checked: the bright soliton sech(x) exp(it/2): its shape, its phase, the conserved norm, and the measured order
4 tests · 71 lines
A nonlinearity against diffusion (nonlinear parabolic)
Steepening or growth on one side, diffusion on the other, and a front whose width and speed are set by the balance.
Burgers' equation, viscous
conservative centered differences in space, fourth-order Runge-Kutta in time
checked: the exact travelling-wave solution, its speed, and the measured order
3 tests · 93 lines
The Fisher-KPP equation
Strang splitting: the logistic reaction solved exactly, diffusion by Crank-Nicolson
checked: the front speed against the KPP value 2, approached from below with Bramson's logarithmic correction
3 tests · 77 lines
The builtins they use
A finite-difference stencil is array arithmetic on shifted copies; an implicit step is a tridiagonal solve; a spectral step is an FFT. These are the builtins the solvers use, with the descriptions from the reference.
roll(a, k) · roll(A, k, axis)- circular shift: element i moves to i + k, wrapping. A periodic stencil in one line: roll(u, 1) - 2.0 * u + roll(u, -1)
tridiag_solve(sub, diag, super, rhs)- the Thomas algorithm: solves a tridiagonal system in O(n), for the diagonally dominant systems that implicit PDE schemes produce. sub[0] and super[n-1] are unused
fft(a)- forward transform, in place, any length
ifft(a)- inverse transform, in place; ifft(fft(x)) returns x to rounding error
slab(A, k)- a COPY of the k-th slab along the leading axis: row k of a matrix, plane k of a rank-3 array. Real or complex
setslab(A, k, S)- overwrite that slab with S, in place
transpose(A)max(a) · min(a)- with ONE array: its largest / smallest entry.
max(abs(u - exact))is the max-norm error of a solve pow(a, s)- every entry of a real array to the same scalar power:
pow(rho, 1.0 / 3.0) sum(a)- the sum of every element of a real array, in NumPy's order in every exact engine: the entries in row-major order, from 0.0, in blocks of 8192, each block by NumPy's pairwise sum (below 8 entries left to right; to 128, eight running sums combined as a tree; above, two halves, the first a multiple of 8 long). Fast mode is not bound by this order
How they are checked
- Against something exact
- A closed-form solution where there is one (a decaying mode, a soliton, d'Alembert, Black–Scholes, the Riemann problem); a manufactured solution where there is not.
- The order is measured, not quoted
- Each solver runs at two or three resolutions and a test asserts the observed convergence rate. That includes the rates below the textbook figure: the first-order shock-capturing scheme for the Euler equations converges at 0.66 in L1, because the contact discontinuity smears, and its test asserts a rate between one half and one.
- Conservation to rounding
- Where the equation conserves something and the scheme is built to, the test tolerance is 1e-12.
- In the interpreter and in C
pde/check.shruns every program through the interpreter and through Viv → C → native, and requires the two outputs to agree line for line.