Stencils and PDE steps

A finite-difference stencil is arithmetic on shifted copies of an array, and an implicit time step is a tridiagonal solve. These are the pieces the PDE solvers are built from. Run a cell to see its output beside the expected one.

A periodic second difference

roll(u, k) moves element i to i + k, wrapping around, so roll(u, 1)[i] is u[i - 1]. The second difference of i squared is exactly 2; the two ends show the wrap.

Viv [1]

u = [0.0, 1.0, 4.0, 9.0, 16.0]
d2 = roll(u, 1.0) - 2.0 * u + roll(u, -1.0)
show d2

roll along one axis of a grid

A third argument picks the axis: 0 shifts rows, 1 shifts along each row. A five-point Laplacian is four of these.

Viv [2]

A = [[1.0, 2.0, 3.0], [4.0, 5.0, 6.0]]
B = roll(A, 1.0, 1.0)
C = roll(A, 1.0, 0.0)
show B[0, 0]
show C[0, 0]

A tridiagonal solve

tridiag_solve(sub, diag, super, rhs) is the Thomas algorithm, O(n). Here the matrix is the 1-D Laplacian (-1, 2, -1), and the right-hand side is chosen so the answer is all ones. sub[0] and super[n-1] are never read.

Viv [3]

n = 4
sub = zeros(n) - 1.0
diag = zeros(n) + 2.0
sup = zeros(n) - 1.0
x = tridiag_solve(sub, diag, sup, [1.0, 0.0, 0.0, 1.0])
show x

One implicit heat step

Backward Euler for u_t = u_xx with dt / h^2 = 1: solve (I - r D2) u_new = u_old. The spike spreads symmetrically, and no time step is too large for it.

Viv [4]

n = 5
r = 1.0
u = [0.0, 0.0, 1.0, 0.0, 0.0]
off = zeros(n) - r
u = tridiag_solve(off, zeros(n) + (1.0 + 2.0 * r), off, u)
show u

The max-norm error

max and min of ONE array reduce it to its largest or smallest entry. With two arguments they stay elementwise.

Viv [5]

x = linspace(0.0, 1.0, 5)
err = max(abs(x * x - x))
show err
show min(x * x - x)

An array to a power

pow(a, s) raises every entry to the same scalar power.

Viv [6]

c = pow([1.0, 8.0, 27.0], 1.0 / 3.0)
show c[1]
show sum(pow([1.0, 2.0, 3.0], 2.0))

Negating an array

Unary minus works on arrays, real or complex.

Viv [7]

a = [1.0, -2.0, 0.5]
b = -a
show b

← Gotchas Small algorithms →