Solving Laplace's equation using finite difference method

Question

Suppose I use a 5-point approximation for the Laplacian $u_{xx} + u_{yy}$ at a point $( x_i , y_j )$ \begin{align*} \frac{-2 u_{ij} + u_{i+1,j} + u_{i-1,j}}{h_x^2} + \frac{-2 u_{ij} + u_{i,j+1} + u_{i,j-1}}{h_y^2} \end{align*} to solve $u_{xx} + u_{yy}=f$ with Dirichlet boundary condition, can I use the Fast Poisson solver? If yes (I think the answer should be yes since the coefficient matrix is of tridiagonal symmetric toeplitz block type), what about imposing a nonzero Dirichlet boundary condition instead to solve $u_{xx} + u_{yy}=0$ using finite difference method? And why the solution vector $u$ cannot have a local maximum at any interior node $(x_{i},y_{j})$? Thanks for any hints.

Answer 1

Consider the discrete Dirichlet problem for Laplace equation: $$ \frac{u_{i-1,j} - 2u_{i,j} + u_{i-1,j}}{h_x^2} + \frac{u_{i,j-1} - 2u_{i,j} + u_{i,j-1}}{h_y^2} = 0, \quad i = 1,\dots,M-1, \quad j=1,\dots,N-1\\ \tag{*} u_{0,j} = g_{0,j}, \quad j=0,\dots,N\\ u_{M,j} = g_{M,j}, \quad j=0,\dots,N\\ u_{i,0} = g_{i,0}, \quad i=1,\dots,M-1\\ u_{i,N} = g_{i,N}, \quad i=1,\dots,M-1 $$ Let's reformulate it as equivalent discrete Possion problem by substituting $v_{ij} = u_{ij} - g_{ij}$. $$ \frac{v_{i-1,j} - 2v_{i,j} + v_{i-1,j}}{h_x^2} + \frac{v_{i,j-1} - 2v_{i,j} + v_{i,j-1}}{h_y^2} = f_{ij}, \quad i = 1,\dots,M-1, \quad j=1,\dots,N-1\\ f_{ij} \equiv -\frac{g_{i-1,j} - 2g_{i,j} + g_{i-1,j}}{h_x^2} - \frac{g_{i,j-1} - 2g_{i,j} + g_{i,j-1}}{h_y^2}\\ v_{0,j} = 0, \quad j=0,\dots,N\\ v_{M,j} = 0, \quad j=0,\dots,N\\ v_{i,0} = 0, \quad i=1,\dots,M-1\\ v_{i,N} = 0, \quad i=1,\dots,M-1 $$ For this to work we need to define $g_{i,j}$ at the internal nodes $i=1,\dots,M-1,j=1,\dots,N-1$.

For direct solvers (direct sparse or FFT methods) the choice for internal $g_{i,j}$ may be arbitrary, for example, $g_{i,j} = 0$. This is equivalent to eliminating boundary unknowns $u_{0,j}$, $u_{M,j}$,$u_{i,0}$, $u_{i,N}$, from the system. After the elimintation the right hand side would contain extra terms containing boundary values of $g_{i,j}$.

For iterative solvers the choice of interval values for $g_{i,j}$ may affect the number of required iterations and it may be wise to take $g_{i,j}$ as a smooth function, for example: $$ g_{i,j} = \tilde c_1 g_{0,j} + \tilde c_2 g_{M,j} + \tilde c_3 g_{i,0} + \tilde c_4 g_{i,N}\\ c_1 = \frac{N h_y}{h_x i}, \quad c_2 = \frac{N h_y}{h_x (M-i)}\\ c_3 = \frac{M h_x}{h_y j}, \quad c_4 = \frac{M h_x}{h_y (N-j)}\\ \tilde c_i = \frac{c_i}{c_1 + c_2 + c_3 + c_4}. $$ This formula is a very rough approximation of Cauchy integral but it gives a pretty smooth $g_{ij}$ if boundary data is also smooth.

The property that $u_{i,j}$ attains its extrema at the boundary is known as discrete maximum principle and does not depend on the way you solve (*), as it is the property of the solution itself. The proof is very simple. Rewrite the discrete Laplace equation as $$ u_{i,j} = \frac{h_y^2 u_{i-1,j} + h_y^2 u_{i+1,j} + h_x^2 u_{i,j+1} + h_x^2 u_{i,j-1}}{2(h_x^2 + h_y^2)}. $$ Assume that $u_{i,j}$ is strictly greater (strictly less) than each of $u_{i-1,j}$, $u_{i+1,j}$, $u_{i,j+1}$, $u_{i,j-1}$. Then the equality cannot hold. So the discrete solution cannot have local maxima or minima inside the domain (that is in points that are surrounded by four neighbours). The only possibility left it that extrema are located at some boundaries.