Laplacian in 2D

Module 4: Space and Pattern Formation

Reading

This laboratory lifts the 1D workflow to a rectangle. The main change is the Laplacian. The Euler loop, the boundary logic, and the validation steps stay the same. If your 1D code was modular, you mostly change the array shape.

Rectangular Domains

You work on \([0, L_x] \times [0, L_y]\). After discretization each field becomes a matrix.

  • The row index stores the \(x\) position, the column index the \(y\) position.
  • The full state has shape (2, num_x, num_y), one matrix for \(u\) and one for \(v\).
num_x = int(length_x / dx)
num_y = int(length_y / dx)
uv = np.ones((2, num_x, num_y))

You still update the whole matrix at once. Shift it up, down, left, and right, and combine the shifted copies.

The 5-Point Stencil

\[ \Delta u = \frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2} \approx \frac{u_{i+1,j} + u_{i-1,j} + u_{i,j+1} + u_{i,j-1} - 4u_{i,j}}{dx^2}. \tag{1}\]

Applying the 1D centered difference in each direction gives this stencil, with error \(\mathcal{O}(dx^2)\) (leveque2007finite?).

Exercise

Implement the stencil with np.roll(). Each call shifts the whole matrix by one grid cell along one direction.

import numpy as np

def laplacian_2d(field, dx):
    # TODO: combine the four nearest neighbours and the centre point with np.roll
    lap =
    return lap / dx**2
def laplacian_2d(field, dx):
    lap = -4 * field
    lap += np.roll(field, shift=1, axis=0)
    lap += np.roll(field, shift=-1, axis=0)
    lap += np.roll(field, shift=1, axis=1)
    lap += np.roll(field, shift=-1, axis=1)
    return lap / dx**2

Boundary Conditions on a Rectangle

The same three options as in Boundary Conditions, now on four edges.

  • Periodic: leave the array unchanged; np.roll already wraps each side onto the opposite one.
  • Neumann: copy the nearest interior row or column onto each edge.
  • Dirichlet: overwrite each edge with prescribed values.
def apply_boundary_conditions(uv, bc="neumann", dirichlet_values=(1.4, 1.96)):
    if bc == "periodic":
        return uv
    if bc == "neumann":
        uv[:, 0, :] = uv[:, 1, :]
        uv[:, -1, :] = uv[:, -2, :]
        uv[:, :, 0] = uv[:, :, 1]
        uv[:, :, -1] = uv[:, :, -2]
        return uv
    if bc == "dirichlet":
        for species, value in enumerate(dirichlet_values):
            uv[species, 0, :] = value
            uv[species, -1, :] = value
            uv[species, :, 0] = value
            uv[species, :, -1] = value
        return uv
    raise ValueError("bc must be 'periodic', 'neumann', or 'dirichlet'")

Validate the Stencil

Test on \(u = \cos(2\pi x/L_x)\cos(4\pi y/L_y)\) with periodic wrapping, where \(\Delta u = -\left[(2\pi/L_x)^2 + (4\pi/L_y)^2\right] u\).

import numpy as np


def laplacian_2d(field, dx):
    lap = -4 * field
    lap += np.roll(field, shift=1, axis=0)
    lap += np.roll(field, shift=-1, axis=0)
    lap += np.roll(field, shift=1, axis=1)
    lap += np.roll(field, shift=-1, axis=1)
    return lap / dx**2


length_x, length_y = 20.0, 50.0
for dx in (1.0, 0.5, 0.25):
    x = np.arange(int(length_x / dx)) * dx
    y = np.arange(int(length_y / dx)) * dx
    xx, yy = np.meshgrid(x, y, indexing="ij")
    u = np.cos(2 * np.pi * xx / length_x) * np.cos(4 * np.pi * yy / length_y)
    exact = -((2 * np.pi / length_x) ** 2 + (4 * np.pi / length_y) ** 2) * u
    print(f"dx = {dx:.2f}   max error = {np.max(np.abs(laplacian_2d(u, dx) - exact)):.2e}")
dx = 1.00   max error = 1.14e-03
dx = 0.50   max error = 2.86e-04
dx = 0.25   max error = 7.15e-05

The error again drops by four for each halving of \(dx\).

Stability in 2D

The 5-point stencil halves the stable time step: \(dt \le dx^2/(4D)\) instead of \(dx^2/(2D)\). With \(d = 30\) and \(dx = 1\) that is \(dt \le 0.0083\), so the 1D choice \(dt = 0.01\) fails in 2D (Time-Step Stability).

Optional extension

A 9-point stencil that also uses the diagonal neighbours gives a more isotropic Laplacian. amlab.pdes_2d.gray_scott.laplacian_9pt implements one. Compare the two stencils on the test function above.

References

Pearson, John E. 1993. “Complex Patterns in a Simple System.” Science 261 (5118): 189–92. https://doi.org/10.1126/science.261.5118.189.