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$ (LeVeque 2007, chap. 1).

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$.

References

LeVeque, Randall J. 2007. Finite Difference Methods for Ordinary and Partial Differential Equations: Steady-State and Time-Dependent Problems. Philadelphia: Society for Industrial; Applied Mathematics. https://doi.org/10.1137/1.9780898717839.

--- title: "Laplacian in 1D" subtitle: "Module 4: Space and Pattern Formation" format: html --- 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, Chapter 1]. ## The Centered Difference $$ \Delta u = \frac{d^2 u}{dx^2} \approx \frac{u(x+dx) - 2u(x) + u(x-dx)}{dx^2}. $$ {#eq-lap1-fd} ::: {.callout-tip} ## 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](time-step-stability.qmd). ::: ## 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](boundary-conditions.qmd)). ::: {.callout-note} ## Exercise Write a helper that computes this second derivative for a 1D field. ```python 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 ``` ::: ::: {.callout-tip collapse="true"} ## Solution ```python 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$. ```{python} 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}") ``` Each halving of $dx$ divides the error by four, as the Taylor argument predicts. ::: {.callout-note} ## 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? ::: ::: {.callout-tip collapse="true"} ## 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$. :::