Discretization of Inhomogeneous Dirichlet Boundary Conditions for 2D Poisson's Equations

Question

Let our problem be $$ \begin{align*} -\Delta u &= f(x,y), \quad (x,y) \in [0,5]\times [0,5]\\ u(x,y) &= g(x,y), \quad (x,y) \in \partial\{[0,5] \times [0,5]\} \end{align*} $$

Suppose you have a $4\times4$ step mesh grid (for a total of $25$ grid points). Begin the indexing at $i,j = (1,1)$ for the bottom left corner; take a lexicographical ordering from bottom up, then progressing from left to right. So in the equation $Au = b$, where $A$ is $25\times25$, b is $25\times1$.

Then for inhomogeneous Dirichlet boundary conditions, what is are the values of $b$ at the corner of the grid and the corner of the interior grid? That is, what values constitute $b_1, b_5, b_{21}, b_{25}$ for the corner of the outer grid, and $b_7, b_9, b_{17}, b_{19}$ for corner points of the interior grid (which has $9$ grid points).

I understand the homogeneous case, but I'm a little confused on the non-homogeneous case.

Answer 1

First, let's enumerate values of $u$ at grid points ($u_{ij}$), so that we have only one index, as stated (lexicographical). This procedure is called flattening of the matrix, i.e. matrix $u_{ij}$ becomes row vector $u_i$. So grid points and function values will be

Simplest laplacian discretization is $$ \Delta u_i = \frac {u_{i-N} + u_{i+N} + u_{i-1} + u_{i+1} - 4u_i}{h^2} = g_i $$ where $N$ is a number of points in each direction, so grid is $N\times N$, and $g_i$ is the function value at grid point where $u_i$ is located. You can check with the picture above to see which points are taken. Obviously, equation makes sense only for the points $u_7$-$u_9$, $u_{12}-u_{14}$, $u_{17}-u_{19}$, since other points are on be boundary and they satisfy BC, not PDE.

So if you flatten your function, for the eligible points $7$-$19$ you get a system of equations. As for the missing points, they satisfy boundary conditions, so $u_i = f_i$, where $f_i$ is the function value where $u_i$ is located. So final system looks as below

Answer 2

I intended to have this as a comment to Kaster's post, but a solution seemed more fit. Note: This is in no way criticizing Kaster's solution or post. In fact, I had the same exact answer as he did, until I decided to investigate further into the wee hours.

The problem with the matrix that is presented in Kaster's solution is that it doesn't handle Dirichlet boundary conditions properly. Notice that you'll in fact double count certain grid nodes, causing the boundary conditions to be Neumann or some sort of mixed boundary condition. Initially, I couldn't get my code$^*$ to run properly -- the $L_2$ error kinked after the first few iterations.

But after fixing the matrix, I was able to get a more accurate result of Matlab: .

Here, the $L_2$ norm displays proper behavior (strictly decreasing); it does slow down once you iterate enough times.

I'm no differential equations expert, but this is what I believe to be correct.

Update

Here's what your sparsity pattern should look like. In your matrix, the positions that are supposed to have values are correct. But you have a few elements that aren't supposed to be there. Think about how the finite difference formulation works and where those elements should be.