Best practice for dealing with Dirichlet boundary conditions in finite-difference schemes: add artificial unknowns?

Question

I know of at least two ways of dealing with Dirichlet boundary conditions in finite-difference schemes (and, to a lesser extent, finite-element schemes). Here I'm thinking of solving Poisson's equation $\Delta u = f$, but the question is more general.

1. Treat the function values at boundary nodes as further unknown variables $u_i$, but add equations of the form $u_i = g(x_i)$ to the resulting system of equations. In total we'll have $n_\text{inner} + n_\text{boundary}$ equations and unknowns, where $n_\text{inner}$ is the number of nodes in the interior of the domain and $n_\text{boundary}$ is the number of nodes on the boundary.

2. Don't include variables for the fixed function values at boundary nodes. In total we'll have $n_\text{inner}$ equations and unknowns.

The two approaches have the following advantages and disadvantages:

• Setting up the equations in the second approach is slightly more complicated, because nodes which are neighbours of boundary nodes have to be treated slightly differently. This seems to be the main selling point of the first approach.

• With the first approach, the resulting system matrix might not be symmetric. (Symmetry is very nice for iterative solvers.)

• The second approach is closer to the theoretical treatment of the undiscretized case: With zero Dirichlet conditions, we look for a solution in the Sobolev space $H_0^1(\Omega)$ instead of $H^1(\Omega)$.

• As soon as more than one differential operator is involved, for instance in the Stokes equation, the first approach gets more messy. One has to introduce zero rows in some of the resulting matrices.

The question: Which of the two approaches is the current recommended best practice? Are there further advantages or disadvantages which I've overlooked?

Answer 1

I think the best approach depends on the solver that you are using. With a direct method, I would incorporate Dirichlet boundary conditions by modifying the right-hand side (option 2). The main advantage, as you have already listed, is that the resulting matrix is symmetric. This not only saves you storage, but also allows for (much) more efficient algorithms.

There are also other advantages. For example, suppose you want to incorporate an irregular boundary in your domain. Option 1 now becomes quite complicated, in particular when the boundary values depend on each other or on the solution -- for example with Neumann or mixed boundary conditions. Option 2 would be the way to go.

With an iterative method, it depends. For example, when using geometric multigrid, it often makes sense to use ghost cells near your domain boundaries so that you can use the same smoother everywhere.