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. All linear, constant-coefficient PDEs with no higher than second derivatives can be classified as either parabolic, hyperbolic, or elliptic. No time variable appears in (13.3.1). Although variable names are arbitrary, elliptic PDEs often do represent systems at steady state.

In order to get a fully specified problem, the Laplace or Poisson equations must be complemented with a boundary condition. Because both $x$ and $y$ are spatial variables, this is our first encounter with a boundary condition that is not imposed simply at a pair of points. 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 use a new matrix operation to solve it.

13.3.2Kronecker product¶

The Kronecker product obeys several natural-looking identities:

13.3.3Poisson 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.7). 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{A}\mathbf{u}=\mathbf{b}$ with the equation

Hence from $\mathbf{A}$ 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 matrix with zeros everywhere except for ones at $(i,i)$ for all $i\in B$. A similar expression can be derived for the modification of $\mathbf{b}$. However, rather than implementing a literal interpretation of (13.3.9), the changes to $\mathbf{A}$ and $\mathbf{b}$ are made more easily through logical indexing. A small example is more illuminating than further description.

13.3.4Implementation¶

Function 13.3.1 is our code to solve the Poisson equation.

We use Demo 13.3.2 as a template: create the linear system, modify it for the boundary conditions, solve it using backslash, and reshape to get a grid function. The matrix is $N=(m+1)(n+1)$ on each side and very sparse, so we take care to use sparse matrices in the code to exploit that structure. There is a small but important change from Demo 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{A}$. This tweak improves the condition number of the final matrix.

13.3.5Accuracy 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 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 sparse Cholesky factorization can be applied to $-\mathbf{A}$, taking $O(n^2N)=O(n^4)$ operations.

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

If we chose a Chebyshev spectral discretization instead, 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. For the specific case of Poisson’s equation on a rectangle, there are specialized fast solution methods that are beyond the scope of this discussion.

13.3.6Exercises¶

1. ✍ Using general $2\times 2$ matrices, verify identities 4 and 6 in Theorem 13.3.1.

2. ✍ Prove that the matrix $\mathbf{A}$ defined in (13.3.7) is symmetric if $\mathbf{D}_{xx}$ and $\mathbf{D}_{yy}$ are symmetric.

3. ⌨ Use Function 13.3.1 to solve the following problems on $[0,1]\times[0,1]$ using $m=n=50$. In each case, make a plot of the solution and a plot of the error.

   (a) $u_{xx}+u_{yy} = 2x^2[x^2(x-1)+(10x-6)(y^2-y)]$, with $u=0$ on the boundary.

   Solution: $u(x,y) = x^4(1-x)y(1-y)$.

   (b) $u_{xx}+u_{yy} = \left(16 x^2 + (1-4 y)^2\right) \sinh (4 x y-x)$, with $u=\sinh(4xy-x)$ on the boundary.

   Solution: $u(x,y) = \sin(4\pi x))$.

   (c) $u_{xx}+u_{yy} = -(20\pi^2) \sin (4\pi x) \cos (2\pi y)$, with $u = \sin (4\pi x) \cos (2\pi y)$ on the boundary.

   Solution: $u(x,y) = \sin (4\pi x) \cos (2\pi y)$.

4. ⌨ For each case in Exercise 3, solve the problem using Function 13.3.1 with $m=n=20,30,40,\ldots,120$. For each numerical solution compute the maximum absolute error on the grid. On a log-log plot, compare the convergence of the error as a function of $n$ to theoretical second-order accuracy.

5. ⌨ Copy Function 13.3.1 to a new function named poischeb, and modify it to use a Chebyshev discretization rather than finite differences. For each item in Exercise 3, solve the problem using poischeb for $m=n=10,15,20,\ldots,40$. For each numerical solution compute the maximum absolute error on the grid. Show the convergence of the error as a function of $n$ on a log-linear plot.

6. ⌨ Sometimes boundary conditions are specified using a piecewise definition, with a different formula for each side of the domain. Use Function 13.3.1 with $m=n=60$ to solve the Laplace equation on $[0,1]^2$ with boundary conditions

   $$u(0,y) = u(1,y) \equiv 0, \quad u(x,0) = \sin(3\pi x), \quad u(x,1) = e^{2x}(x-x^2).$$

   (13.3.10)

   Make a surface plot of your numerical solution.