# Dirichlet Boundary Condition finite difference method using sparse-matrix $Ax = b$ system

I am trying to solve the boundary value problem for heat equation:

$$u_{xx} + u_{yy} = f(x,y)$$

where the solution $$u(x,y) \in [0,1] \times [0,1]$$ and the Dirichlet boundary condition $$u(x,y) = u_0$$.

So I discretized the problem with $$N$$ interior points with step size $$h = 1/(N+1)$$

$$w_{i,j-1} + w_{i+1, j} + w_{i-1,j} - 4 w_{i,j} + w_{i, j+1} = h^{2} f(x_{i}, y_{j})$$

where $$i,j \in \{1,2,...,N\}$$ and the boundary conditions are:

$$w_{0,j} = w_{N+1, j} = u_{0} \quad\quad\quad w_{i,0} = w_{i,N+1} = u_{0}$$

where $$i,j$$ is from $$1$$ to $$N$$

I reduce this problem into sparse-matrix system. I know how to make the sparse matrix $$A$$ for this problem by using Python. The only issue is the RHS which is matrix $$b$$. For the case of $$u_0 = 0$$, I trivially evaluate the grid points at function $$f(x,y)$$. For the case $$u_0 \neq 0$$, I do not know to retrieve matrix $$b$$. I tried to write down with the case $$N = 3$$ to see the general pattern as below (I put in column-order):

I hope anyone could help me understand how to retrieve that matrix $$b$$ in general case for Dirichlet boundary condition??. The code I wrote for this problem is in Python, but I want to understand the generality first.

• It may be useful to think about it in 1D first. The matrix equation at a boundary grid point $(i,j)$ should express the relation $w_{i,j}=u_0$, so the corresponding entry in matrix $A$ should be 1, and in vector $b$ is should be $u_0$. – Maxim Umansky Oct 3 '20 at 2:20
• But then it would be totally different from this problem. Isn't it? The discretization in 1D problem is taking 3 stencils point. Meanwhile in this, I think it is 5 stencils (grid points) – Dong Le Oct 3 '20 at 2:25
• I see what you are saying. If $A = I$ which is the identity matrix, then you are left with $\vec{w} = \vec{b}$ where $\vec{b} = h^{2} \vec{f} - \vec{U_{0}}$. So unless all $f_{i,j} = 0$ then you have the relation $w_{i,j} = - U_{i,j}$ – Dong Le Oct 3 '20 at 2:32
• On the boundary $f$ should not be present, only the boundary condition relation $w_{i,j}=u_0$ – Maxim Umansky Oct 3 '20 at 4:05

A general way would be to include the boundary nodes in the definition of $$A$$ (which will give you a matrix with more columns than rows) and derive $$b$$ as the contribution of the Dirichlet nodes. This way, other linear terms like convection are readily included.

Assume that the matrix $$A$$ looks like

$$A = [A_I | A_\Gamma ]$$

where $$A_I$$ is the (square) operator on the inner (as you have it in your question) and where $$A_\Gamma$$ are the columns that belong to the boundary nodes. It is not a problem to resort the columns like this, but in python I would rather work with the indices and slices.

Then the $$b$$ is extracted as

$$b = - A [0 \dotsm 0|u_0 \dotsm u_0]^T = -A_\Gamma[u_0 \dotsm u_0]^T.$$

In fact, your solution reads $$w=[w_1 \dotsm w_{N^2} | u_0 \dotsm u_0]$$ and will solve $$Aw=f \quad \text{ or } \quad A_Iw_I + A_\Gamma w_\Gamma = f \quad \text{ or }\quad A_Iw_I = -A_\Gamma w_\Gamma + f \quad \text{or}\quad A_Iw_I = b + f ,$$ where $$w_\Gamma$$ is your solution at the boundary nodes. We have once explained this approach for Finite Elements discretizations in a preprint^1 but the principle is the same for Finite Differences.