Laplace and Poisson equations

Consider the heat equation $u_t=u_{xx}+u_{yy}$. After a long time, the distribution of temperature will stop changing. This steady-state solution must satisfy the PDE $u_{xx}+u_{yy}=0$, which is our third and final canonical PDE.

The Laplace/Poisson equation is the archetype of an elliptic PDE, which is our third and final category for linear PDEs of second order; we have already encountered parabolic (i.e., diffusion) and hyperbolic (advection and wave) examples.

In order to get a fully specified problem, the Laplace or Poisson equation must be complemented with a boundary condition. We consider only the Dirichlet condition $u(x,y)=g(x,y)$ around the entire boundary.

13.3.1Sylvester equation¶

With the unknown solution represented by its values $\mathbf{U}=\mtx(u)$ on a rectangular grid, and second-derivative finite-difference or spectral differentiation matrices $\mathbf{D}_{xx}$ and $\mathbf{D}_{yy}$, the Poisson equation (13.3.1) becomes the discrete equation

where $\mathbf{F}=\mtx(f)$. Equation (13.3.2), with an unknown matrix $\mathbf{U}$ multiplied on the left and right in different terms, is known as a Sylvester equation. We will introduce a new matrix operation to solve it.

The Kronecker product obeys several natural-looking identities:

13.3.2Poisson as a linear system¶

We can use (13.3.4) to express the Sylvester form (13.3.2) of the discrete Poisson equation as an ordinary linear system. First, we pad (13.3.2) with identity matrices,

where $\mathbf{I}_{x}$ and $\mathbf{I}_{y}$ are the $(m+1)\times (m+1)$ and $(n+1)\times (n+1)$ identities, respectively. Upon taking the vec of both sides and applying (13.3.4), we obtain

This is in the form of a standard linear system in $(m+1)(n+1)$ variables.

Boundary conditions of the PDE must yet be applied to modify (13.3.6). As has been our practice for one-dimensional boundary-value problems, we replace the collocation equation for the PDE at each boundary point with an equation that assigns that boundary point its prescribed value. The details are a bit harder to express algebraically in the two-dimensional geometry, though, than in the 1D case.

Say that $N=(m+1)(n+1)$ is the number of entries in the unknown $\mathbf{U}$, and let $B$ be a subset of $\{1,\dots,N\}$ such that $i\in B$ if and only if $(x_i,y_i)$ is on the boundary. Then for each $i\in B$, we want to replace row $i$ of the system $\mathbf{L}\mathbf{u}=\mathbf{f}$ with the equation

Hence, from $\mathbf{L}$ we should subtract away its $i$th row and add $\mathbf{e}_i^T$ back in. When this is done for all $i\in B$, the result is the replacement of the relevant rows of $\mathbf{A}$ with relevant rows of the identity:

where $\mathbf{I}_B$ is a square matrix with zeros everywhere except for ones at $(i,i)$ for all $i\in B$. A similar expression can be written for the modification of $\mathbf{f}$,

where $\mathbf{g} = \operatorname{vec}(\operatorname{mtx}(g))$ is the vector of boundary function evaluations. The result is a modified linear system $\mathbf{A}\mathbf{u}=\mathbf{b}$ that has incorporated the desired boundary conditions.

While (13.3.8) and (13.3.9) are algebraically correct, a literal implementation of them is not the most efficient route. Instead, the relevant changes are made more easily through logical indexing into rows and columns. A small example is more illuminating than further description.

13.3.3Implementation¶

Function 13.3.1 is our code to solve the Poisson equation.

We use Example 13.3.2: the boundary conditions are rescaled to read $\sigma u(x,y)=\sigma g(x,y)$, where σ is the largest element of a row of $\mathbf{L}$. This tweak improves the condition number of the final matrix.

13.3.4Accuracy and efficiency¶

In Function 13.3.1 we used second-order finite differences for the discretization. Let’s simplify the discussion by assuming that $m=n$. The truncation error in the solution can be expected to be $O(n^{-2})$.

The matrix $\mathbf{A}$ has size $N=O(n^2)$. The upper and lower bandwidths of the matrix $\mathbf{A}$ are both $O(n)$. It can be shown that the matrix is symmetric and negative definite, so banded Cholesky factorization can be applied, taking $O(n^2N)=O(n^4)$ operations.

As $n$ increases, the truncation error decreases while the operation count increases. In fact, the growth in operations is faster than the corresponding decrease in error, making it costly to get high accuracy. Say we have a fixed running time $T$ at our disposal. Then $n=O(T^{1/4})$, and the convergence of error is $O(1/\sqrt{T})$. For instance, reduction of the error by a factor of 10 requires 100 times the computational effort.

If we chose a Chebyshev spectral discretization instead of finite differences, the calculus changes. Provided that the solution is smooth, we expect convergence at a rate $K^{-n}$ for some $K>1$. However, the system matrix is no longer sparse nor symmetric, and solution by LU factorization takes $O(N^3)=O(n^6)$ flops. Hence, as a function of running time $T$, we expect a convergence rate on the order of $K^{-T^{1/6}}$. It’s not simple to compare this directly to the finite-difference case. In the long run, the spectral method will be much more efficient, but if the accuracy requirement is undemanding, second-order finite differences may be faster.

13.3.5Exercises¶