The Discretized Laplace Operator on Hyperrectangles with Zero Dirichlet Boundary Conditions

In this blog post, I present stiffness and mass matrix as well as eigenvalues and eigenvectors of the Laplace operator (Laplacian) on domains $(0, \ell)$, $(0, \ell_1) \times (0, \ell_2)$, and so on (hyperrectangles) with zero Dirichlet boundary conditions discretized with the finite difference method (FDM) and the finite element method (FEM) on equidistant grids. For the FDM discretization, we use the central differences scheme with the standard five-point stencil in 2D. For the FEM, the ansatz functions are the hat functions. The matrices, standard eigenvalue problems $A v = \sigma v$, and generalized eigenvalue problems $K w = \tau M w$ arising from the discretization lend themselves for test problems in numerical linear algebra because they are well-conditioned, not diagonal, and the matrix dimension can be increased arbitrarily.

Python code generating the matrices and their eigenpairs can be found in my git repository discrete-laplacian.

Introduction

Let $d \in \mathbb{N} \setminus \{0\}$, let $u: \mathbb{R}^d \rightarrow \mathbb{R}$. The differential operator $\Delta$ is called Laplacian and it is the sum of the second derivatives of a function:

It is a well-researched differential operator in mathematics and for many domains, the exact eigenvalues $\lambda$ and eigenfunctions $u$ are known such that

In this blog post, I discuss the solutions of $-\Delta u = \lambda u$ on hyperrectangles $\Omega$ with the Dirichlet boundary condition

The combination of Laplace operator and the finite difference method (FDM) can be very well used in an introductory course on the numerical treatment of partial differential equations (PDEs) for the illustration of concepts such as discretization of PDEs, discretization error, the growth of the matrix condition number with finer grids, and sparse solvers.

In numerical linear algebra, the Laplace operator is appealing because the FDM discretization of the operator on a one-dimensional domain yields a standard eigenvalue problem $A v = \sigma v$ with a sparse, real symmetric positive-definite, tridiagonal Toeplitz matrix $A$ and known eigenpairs. For domains in higher dimensions, the matrices can be constructed with the aid of the Kronecker sum and consequently, the eigenpairs can be calculated in this case, too. With this knowledge, the Laplace operator makes for a good, easy test problem in numerical mathematics because we can distinguish between discretization errors and algebraic solution errors. Naturally, the Laplace operator can also be discretized with the finite element method (FEM) yielding a generalized eigenvalue problem $Kw = \tau Mw$ with sparse, real symmetric positive-definite matrices. Here, the eigenpairs are known, too, and furthermore, the eigenvectors are exact in the grid points (so-called superconvergence) providing us with another well-conditioned test problem.

I will present the solutions for the continuous case before I briefly introduce Kronecker products and Kronecker sums since these linear algebra operations are used to construct the matrices corresponding to higher dimensional domains. Finally, I discuss domain discretization and notation before giving closed expressions for the matrices and their eigenpairs created by FDM and FEM. In the end, there is an example demonstrating the use of the quantities $\ell_d$, $n_d$, and $h_d$.

Continuous Case

In 1D, the exact eigenpairs $(\lambda_i, u_i(x))$ of the Laplace operator on the domain $(0, \ell)$ are

In 2D, the exact eigenpairs $(\lambda_{ij}, u_{ij}(x))$ of the Laplace operator on the domain $(0, \ell_1) \times (0, \ell_2)$ are

On a $d$-dimensional hyperrectangle, the eigenpairs are

These solutions can be found, e.g., in the book Methods of Mathematical Physics, Vol. I, Chapter VI, §4.1.

Kronecker Products

The matrices of the discretized Laplacian on higher dimensional domains can be constructed with Kronecker products. Given two matrices $A = [a_{ij}] \in \mathbb{R}^{m,n}$, $B \in \mathbb{R}^{k,\ell}$, the Kronecker product $C = A \otimes B$ is defined as follows:

If $A$ and $B$ are square, then $C$ is square as well; let $(\lambda_i, v_i)$ be the eigenpairs of $A$ and let $(\mu_j, w_j)$ be the eigenpairs of $B$. Then the eigenpairs of $C$ are

If $A$ and $B$ are real symmetric, then $C$ is also real symmetric. For square $A$ and $B$, the Kronecker sum of $A$ and $B$ is defined as

where $I_n$ is the $n \times n$ identity matrix. The eigenvalues of $D$ are

The Kronecker sum occurs during the construction of the 2D FDM matrix. See, e.g., Matrix Analysis for Scientists and Engineers by Alan J. Laub, Chapter 13, for more information on these operations.

Domain Discretization

For the 1D case along the $d$-th axis, we use $n_d+2$ points uniformly distributed over $(0, \ell_d)$, such that the step size is $h_d = \frac{\ell_d}{n_d+1}$. We use lexicographical ordering of points, e.g., in 2D let $u: \mathbb{R}^2 \rightarrow \mathbb{R}$ be an eigenfunction of the continuous problem, let $v \in \mathbb{R}^{n_1 n_2}$ be an eigenvector of the algebraic eigenvalue problem. Then $v_1$ is an approximation to $u(x_{1,1})$, $v_2$ approximates $u(x_{2,1})$, $v_{n_1+1}$ approximates $u(x_{1,2})$, and so on.

To obtain accurate approximations to the solutions of the continuous eigenvalue problem, the distance between adjacent grid points should always be constant. Thus, if the length of the 1D domain is doubled, the number of grid points should be doubled, too.

Notation

In the following sections, we need to deal with matrices arising from discretizations of the Laplace operator on 1D domains with different step sizes and their eigenpairs as well as discretizations of the Laplace operator on higher dimensional domains and their eigenpairs. $I$ denotes the identity matrix; its dimension can be gathered from the context.

Matrices $A$ denote the Laplacian discretized with the FDM, matrices $A^{(d)} \in \mathbb{R}^{n_d,n_d}$ denote the discrete Laplacian on 1D domains along the $d$-th axis, and matrices $A_d$ denote the the discrete Laplacian on a $d$-dimensional domain. The eigenpairs of $A^{(d)}$ are signified by $(\sigma_i^{(d)}, v_i^{(d)})$, $i = 1, 2, \dotsc, n_d$, $(\sigma_{i_1 i_2 \dotsb i_d}, v_{i_1 i_2 \dotsb i_d})$ are the eigenpairs of $A_d$.

Similarly, $K$ and $M$ are the stiffness and mass matrix, respectively, of the Laplace operator discretized with the finite element method. $K^{(d)}, M^{(d)} \in \mathbb{R}^{n_d,n_d}$ denote the discrete Laplacian on a one-dimensional domain along the $d$-th axis, and $K_d$, $M_d$ are the stiffness and mass matrix for the discrete Laplacian on a $d$-dimensional hyperrectangle. We speak of the solutions of $Kw = \tau M w$ or of eigenpairs of the matrix pencil $(K, M)$. The eigenpairs of $(K^{(d)}, M^{(d)})$ are denoted by $(\tau_i^{(d)}, w_i^{(d)})$, whereas $(\tau_{i_1 i_2 \dotsb i_d}, w_{i_1 i_2 \dotsb i_d})$ signify eigenpairs of $(K_d, M_d)$.

Finite Difference Method

In this section, we will construct the matrices of the discretized Laplace operator on a $d+1$-dimensional domain with the aid of the matrices for the $d$-dimensional and the one-dimensional case.

Let $e_i$ be the $i$-th unit vector. Let us consider the discretization along the $d$-th axis first. Let

be the forward difference quotient, let

be the backward difference quotient. Then

Consequently, the discrete Laplacian has the stencil

meaning $A^{(d)}$ is a real symmetric tridiagonal $n_d \times n_d$ matrix with the value 2 on the diagonal and -1 on the first sub- and superdiagonal. The eigenvalues of $A^{(d)}$ are

the entries of the eigenvector $v_i^{(d)}$ are

In 2D with lexicographic ordering of the variables, we have

and this matrix can be constructed as follows:

The eigenpairs can be derived directly from the properties of the Kronecker sum: the eigenvalues are

and the eigenvectors are

where $i = 1,2, \dotsc, n_1$, $j=1, 2, \dotsc, n_2$.

In higher dimensions, it holds that

where

are the eigenvalues and

are the eigenvectors, $i_{d+1} = 1, 2, \dotsc, n_{d+1}$.

Finite Element Method

Similarly to the finite differences method, we can construct the matrices of the discretized equation $-\Delta u = \lambda u$ recursively. We use hat functions as ansatz functions throughout this section.

The 1D stencil is

for the stiffness matrix and

for the mass matrix. The generalized eigenvalues are

and the $k$-th entry of the eigenvector $w_i^{(d)}$ corresponding to $\tau_i^{(d)}$ is

Observe that $w_i^{(d)}$ is not only an eigenvector of the matrix pencil $(K^{(d)}, M^{(d)})$ but also of the matrices $K^{(d)}$ and $M^{(d)}$ themselves. Thus, the eigenvalues of $K^{(d)}$ are

whereas the eigenvalues of $M^{(d)}$ are

In 2D, the mass matrix stencil is

Judging by the coefficient $\frac{1}{36} h_1 h_2$ and the diagonal blocks, it must hold that

The stiffness matrix stencil is

Seeing the factors $h_1$ and $h_2$, the stiffness matrix cannot be the Kronecker product or the Kronecker sum of 1D stiffness matrices. However, observe that the 1D mass matrices $M^{(1)}$ and $M^{(2)}$ have coefficients $h_1$ and $h_2$, respectively, and indeed, the 2D stiffness matrix can be constructed with the aid of the 1D mass matrices:

Based on the properties of the Kronecker product and using the fact that mass and stiffness matrix have the same set of eigenvalues, $(K_2, M_2)$ has the eigenvalues

and eigenvectors

where $i = 1, 2, \dotsc, n_1$, $j = 1, 2, \dotsc, n_2$.

For the Laplacian on hyperrectangles in $d+1$-dimensional space, the stiffness matrix can be constructed by

the mass matrix can be constructed using

such that $(K_{d+1}, M_{d+1})$ has eigenvalues

and eigenvectors

where $i_{d+1} = 1, 2, \dotsc, n_{d+1}$.

Example

Now I will present the matrices of the discrete Laplacian on the domain $(0, 6) \times (0, 5)$. The figure below shows the domain and the grid used for discretization: there are five interior grid points along the first axis, three interior points along the second axis, and the interior grid points are highlighted as black circles. Thus, $\ell_1 = 6$, $\ell_2 = 5$, $n_1 = 5$, $n_2 = 3$, the step size $h_1 = 1$, and $h_2 = 1.25$. In 2D, an eigenvector $v$ of the algebraic eigenvalue problem will possess entries $v_i$, $i = 1, 2, \dotsc, 15$, such that $v_1 \approx u(x_{1,1})$, $v_2 \approx u(x_{2,1})$, $\dotsc$, $v_5 \approx u(x_{5,1})$, $v_6 \approx u(x_{1,2})$, $\dotsc$, $v_{15} \approx u(x_{5,3})$ because of the lexicographic ordering. An equispaced grid on a 2D domain with five interior grid points along the first dimension and three interior grid points along the second dimension

In the first dimension, the Laplacian discretized with the FDM has the following matrix:

Along the second axis, we have

With the FEM, the stiffness matrix is

and the mass matrix is

along the first axis. The matrices corresponding to the second dimension are

and