TransWikia.com

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

Computational Science Asked on April 14, 2021

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} quadquadquad 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.

One Answer

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.

Answered by Jan on April 14, 2021

Add your own answers!

Ask a Question

Get help from others!

© 2024 TransWikia.com. All rights reserved. Sites we Love: PCI Database, UKBizDB, Menu Kuliner, Sharing RPP