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):

enter image description here

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.

  • $\begingroup$ 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$. $\endgroup$ – Maxim Umansky Oct 3 '20 at 2:20
  • $\begingroup$ 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) $\endgroup$ – Dong Le Oct 3 '20 at 2:25
  • $\begingroup$ 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}$ $\endgroup$ – Dong Le Oct 3 '20 at 2:32
  • 1
    $\begingroup$ On the boundary $f$ should not be present, only the boundary condition relation $w_{i,j}=u_0$ $\endgroup$ – 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.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.