Laplacian in 1D

Module 4: Space and Pattern Formation

The first implementation step is the spatial derivative. In one dimension the Laplacian is the second derivative, and you approximate it with a centered finite difference on a grid with spacing \(dx\) (leveque2007finite?).

The Centered Difference

\[ \Delta u = \frac{d^2 u}{dx^2} \approx \frac{u(x+dx) - 2u(x) + u(x-dx)}{dx^2}. \tag{1}\]

Where does this formula come from?

Expand \(u\) around \(x\) with Taylor’s theorem:

\[ \begin{aligned} u(x+dx) &= u(x) + u'(x)\,dx + \tfrac{1}{2} u''(x)\,dx^2 + \tfrac{1}{6} u'''(x)\,dx^3 + \mathcal{O}(dx^4), \\ u(x-dx) &= u(x) - u'(x)\,dx + \tfrac{1}{2} u''(x)\,dx^2 - \tfrac{1}{6} u'''(x)\,dx^3 + \mathcal{O}(dx^4). \end{aligned} \]

Adding the two lines cancels the odd terms:

\[ \frac{u(x+dx) - 2u(x) + u(x-dx)}{dx^2} = u''(x) + \mathcal{O}(dx^2). \]

The error is second order in \(dx\): halving \(dx\) divides it by about four. You will test this on Time-Step Stability.

Vectorize with np.roll

On the grid \(x_i = i\,dx\), the formula uses the neighbours \(u_{i-1}\) and \(u_{i+1}\). np.roll(values, 1) shifts the array one place to the right, so position \(i\) holds \(u_{i-1}\). Rolling wraps the ends around, which is exactly a periodic boundary. For other boundaries you overwrite the edge values after each step (Boundary Conditions).

Exercise

Write a helper that computes this second derivative for a 1D field.

import numpy as np

def laplacian_1d(values, dx):
    # TODO: use the left and right neighbors of each grid point
    lap =
    return lap / dx**2

Solution

def laplacian_1d(values, dx):
    lap = -2 * values
    lap += np.roll(values, shift=1)
    lap += np.roll(values, shift=-1)
    return lap / dx**2

Validate on a Known Function

Test the helper on \(u(x) = \cos(2\pi m x / L)\), a periodic function with \(u'' = -(2\pi m/L)^2 u\).

import numpy as np


def laplacian_1d(values, dx):
    lap = -2 * values
    lap += np.roll(values, shift=1)
    lap += np.roll(values, shift=-1)
    return lap / dx**2


length, m = 40.0, 3
for num_points in (40, 80, 160):
    dx = length / num_points
    x = np.arange(num_points) * dx
    u = np.cos(2 * np.pi * m * x / length)
    exact = -((2 * np.pi * m / length) ** 2) * u
    error = np.max(np.abs(laplacian_1d(u, dx) - exact))
    print(f"dx = {dx:.3f}   max error = {error:.2e}")

dx = 1.000   max error = 4.08e-03
dx = 0.500   max error = 1.03e-03
dx = 0.250   max error = 2.57e-04

Each halving of \(dx\) divides the error by four, as the Taylor argument predicts.

Exercise

The course package applies the stencil to both species at once with np.roll(uv, 1, axis=1), where uv has shape (2, N). Why is axis=1 needed?

Solution

Axis 0 indexes the species \(u\) and \(v\), axis 1 indexes space. Without axis=1, NumPy flattens the array, so the last point of \(u\) would become the left neighbour of the first point of \(v\).