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

Question

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.

Answer 1

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.