Consider $\mathcal{L}u(x) = -\Delta u(x,y) = b(x,y)$ on the domain $\Omega = (0,1)^2$ with boundary conditions $u(x,y) = g(x,y)$.
Exercise: Set-up the discrete problem (matrix and right-hand side) for the standard five-point discretization of $\mathcal{L}$ and $h = \frac{1}{3}$
a) without elimination of boundary conditions.
b) with elimination of boundary conditions.
using a lexicographic ordering of grid points in both.
Question: how do I handle/solve this?
Thanks in advance!
On
If $(O,x,y)$ is a Cartesian coordinate system of the plane, then the Laplace operator is given by $-\mathcal{L}u = u_{xx} + u_{yy}$. Using centered finite differences of the second derivatives, we have $$ -\mathcal{L}_hu = \frac{u_{i+1,j} - 2u_{i,j} + u_{i-1,j}}{h^2} + \frac{u_{i,j+1} - 2u_{i,j} + u_{i,j-1}}{h^2} \, , $$ where $u_{i,j} \simeq u(ih,jh)$. Therefore, if $b_{i,j} = b(ih,jh)$, the equation $-\mathcal{L}_hu = b$ becomes $$ -\frac{u_{i-1,j} + u_{i,j-1} - 4u_{i,j} + u_{i,j+1} + u_{i+1,j}}{h^2} = b_{i,j} \, $$ in the interior of the domain where the indices are $0<i,j<1/h$. At the boundaries, we have $u_{i,j} = g_{i,j} = g(ih,jh)$, with the index $i$ or $j$ in $\lbrace 0,1/h\rbrace$. This system of linear equations can be written in the form ${\bf A}{\bf X} = {\bf B}$, where $\bf X$ is the vector of unknowns $u_{i,j}$ sorted in lexicographical order.
a) Without elimination of the boundary conditions, the vector of unknowns ${\bf X} = (u_{i,j})^\top_{0\leq i,j\leq 1/h}$ is in $\Bbb{R}^{(1/h+1)^2}$. We write $(1/h-1)^2$ interior equations and $4/h$ boundary conditions using a square matrix $\bf A$ of $\mathcal{M}_{(1/h+1)^2}(\Bbb{R})$ and a vector $\bf B$ of $\Bbb{R}^{(1/h+1)^2}$. These matrices have the following structure: $$ {\bf A} = \left( \begin{array}{cccc} {\bf I} & & & & \\ {\bf 0} & {\bf E}^{(1,1)} & \dots & {\bf E}^{(1/h-1,1)} & {\bf 0} \\ \hline & & {\bf L} & & \\ \hline {\bf 0} & {\bf E}^{(1,1/h+1)} & \dots & {\bf E}^{(1/h-1,1/h+1)} & {\bf 0} \\ & & & & {\bf I} \end{array}\right) , \quad {\bf B} = \left(\begin{array}{c} (g_{0,j})_{0\leq j\leq 1/h}^\top\\ (g_{i,0})_{0<i<1/h}^\top\\ \hline (b_{i,0<j<1/h})_{0<i<1/h}^\top\\ \hline (g_{i,1/h})_{0<i<1/h}^\top\\ (g_{1/h,j})_{0\leq j\leq 1/h}^\top\end{array}\right) , $$ where $\bf I$ and $\bf 0$ are respectively the identity and zero matrices of $\mathcal{M}_{1/h+1}(\Bbb{R})$. The coordinates of the matrix ${\bf E}^{(a,b)}$ with $1/h-1$ rows ans $1/h+1$ columns are equal to one at the position $(a,b)$ and zero elsewhere. The matrix $\bf L$ is deduced from $\mathcal{L}_h$. In the case $h = 1/2$, we have $$ {\bf A} = \left( \begin{array}{c} 1 & \\ & 1 & \\ & & 1 \\ & & & 1 \\ \hline & -1/h^2 & & -1/h^2 & 4/h^2 & -1/h^2 & & -1/h^2 & \\ \hline & & & & & 1 & & \\ & & & & & & 1 & \\ & & & & & & & 1 & \\ & & & & & & & & 1 \end{array}\right) . $$
b) With elimination of the boundary conditions, the discrete equations satisfied at the boundaries are injected in $\mathcal{L}_h$ at the grid nodes $(i,j)$, where $i$ or $j$ is in $\lbrace 1,1/h-1\rbrace$. Therefore, the vector of unknowns reduces to ${\bf X} = (u_{i,j})^\top_{0< i,j< 1/h}$ which belongs to $\Bbb{R}^{(1/h-1)^2}$, and the matrices $\bf A$ and $\bf B$ in the system ${\bf A}{\bf X} = {\bf B}$ are modified. In the case $h = 1/2$, we are left with $$ {\bf A} = 4/h^2 \, ,\quad {\bf B} = b_{1,1} + (g_{0,1} + g_{1,0} + g_{1,2} + g_{2,1}) / h^2 , $$ so that the solution of the linear system is $u_{1,1} = \frac{1}{16} b_{1,1}+ \frac{1}{4}(g_{0,1} + g_{1,0} + g_{1,2} + g_{2,1})$.
Note: the solution is detailed in the case $h=1/2$ for simplicity purposes, but the same method applies for $h = 1/3$.
Imagine a square grid with inner grid distance $1/3$. Since we are working on $(0,1)^2$ we this grid will consist of $4 \times 4$ grid points.
We start numbering them row by row: the left bottom corner corresponding to $u(0,0)$ is first, the right bottom corner $u(1,0)$ will be fourth and so on. See the image below.
Now, lets figure out what must hold for the different grid points. The outer "ring" of points corresponds to the BC and there we simply want $u(x,y) = g(x,y)$. In the discrete scheme (which will consist of a $4^2 \times 4^2$ matrix A, $4^2 \times 1$ vector $u$ we want to solve for and RHS vector $b$) we can set a $1$ on the main diagonal in the row corresponding to the point number $k$ and the we set $k$'th value of $b$ equal to $g(x,y)$.
For example: 5'th row of matrix $A$ will read $[0, 0, 0, 0, 1, 0, .... 0]$ and 5'th entry of vector $b$ would be $g(0, 1/3)$.
For the interior points we must have: $-\nabla^2 u = b$. We approximate $\nabla^2 u$ with the $5$ point stencil. It states:
$$\nabla^2 u_{i,j} \approx \frac{u_{i-1,j} + u_{i+1,j} + u_{i,j-1} + u_{i,j+1} - 4 u_{i,j}}{h^2}$$
Where $i,j$ correspond to rows and columns of the grid. How to translate $i,j$ into our lexicographical ordering? Well, neighbors to the left and right ($u_{i-1,j}, u_{i+1,j}$ are easy enough - if we are on node $k$ then $k-1, k+1$ are nodes you are looking for. To get a step up or down we need to go through a whole row - so in our case $4$. Thus, we are looking at $k - 4, k+4$.
The row of $A$ corresponding to a interior point should therefore read: $\frac{1}{h^2} [0, ..., 1, 0, 0, 1, -4, 1, 0, 0, -1, 0.... 0]$ and the RHS side vector will read $b_{i,j}$.
We are now only left with figuring out which nodes are interior, which ones are boundaries. This is lucky easy to do: first $4$ nodes ($\frac{1}{h} +1$) are boundaries, then you get one boundary, $2$ ($\frac{1}{h} -1$) interior and again a boundary. This is repeated until we finish with 4 boundary nodes.
Thus, we will have: $5-2-2-2-5$ pattern of boundary/interior nodes. We can now write down the set of equations as shown below.
You should note something here: the calculations done on the boundary points are trivial. Furthermore, if $h$ becomes very small you the entries of your matrix will get small and this can lead to issues. We will now try to exclude the boundary points.
As mentioned, we have $(\frac{1}{h}-1)^2 = 4$ interior points. Thus, our matrix will only be $4\times4$. Lets for example write out the 6'th equation of our system. $$\frac{1}{h^2}(u(\frac{1}{3},0) + u(0,\frac{1}{3}) - 4u(\frac{1}{3},\frac{1}{3}) + u(\frac{1}{3}, \frac{2}{3}) + u(\frac{2}{3},\frac{2}{3})) = b(\frac{1}{3},\frac{1}{3})$$
The thing is, since our bc are trivial we actually know the first two terms! They must equal, $\frac{1}{h^2} (g(\frac{1}{3},0) + g(0,\frac{1}{3}))$. So, we set them to the RHS and are left with: $$\frac{1}{h^2}(-4u(\frac{1}{3},\frac{1}{3}) + u(\frac{1}{3}, \frac{2}{3}) + u(\frac{2}{3},\frac{2}{3})) = b(\frac{1}{3},\frac{1}{3})-\frac{1}{h^2} (g(\frac{1}{3},0) + g(0,\frac{1}{3}))$$
We do it with every interior point that has a boundary point as an neighbor. (in this example - that would be every interior point).
We are left with: