Điểm biên trong Finite Difference cho phương trình 2-D Numerical PDE

[PhuongDong0701 (Phương Đông)]

• Chào anh chị, phần em đang implement cách tính finite difference cho phương trình PDE bậc 2. Cụ thể em có phương trình heat:

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

g(x,y) = 5000\exp(-(x-1/4)^2 + (y-3/4)^2)) với điểm biên là u(x,y) = 300 trên \Delta(D), D(\Omega) = [0,1]\times[0,1] là miền xác định solution của em. Phần là em đã giải ra được cái ma trận Laplacian (discretization):

w_{i+1,j} - 4w_{ij} + w_{i-1,j} + w_{i,j+1} + w_{i,j-1} = h^{2} g(x_i, y_j)

với w_{0,j} = w_{N+1, j} = 1 (j = 1,...,N) ; w_{i,0} = w_{i,N+1} = 1 (i = 1,...,N)

Em có cái sparse matrix cho Laplacian với cái h^{2}g(x_i,y_j). Nhưng em ko biết làm sao generate cái ma trận theo điểm biên

"""The function constructs a sparse matrix A 
    containing the 2D Laplacian matrix
    The dimension of matrix A is n = N^2
    """

def build_laplace_2D(N): 
    n = N*N; 
    
    main_diags = np.ones(n)*-4.0 ## The main -4's diagnonals. 
    
    side_diags = np.ones(n-1) ## Side 1 diagnonals. 
    
    side_diags[np.arange(1,n)%4==0] = 0
    
    up_down_diags = np.ones(n-3)
    
    diagonals = [main_diags,side_diags,side_diags,up_down_diags,up_down_diags]; 
    
    ### Create laplacian matrix: 

    laplacian = sparse.diags(diagonals, [0, -1, 1,N,-N], format="csr")
    Output = laplacian.toarray(); ## convert to np array. 
    
    return Output


N =  64; #  interior point

n = N * N
h = 1 / (N + 1) # step size
h2 = 1/((N+1)**2) # 


delta_x = h*np.arange(1, N+1)
delta_y = h*np.arange(1, N+1)[:, None]

### Calculate f in row lexicographic order #### 
calcul_f2= h2*G(delta_x, delta_y).ravel(order='F')
feval2 = np.reshape(calcul_f2, (n,1), order  =  "F") 

### Generate the  2d Laplacian matrix: 
A = build_laplace_2D(N) 
print(A.shape)
Ainv = np.linalg.inv(A)

u_numerical_lexi2  = np.matmul(Ainv,feval2, order = "F")
u_numerical2 = np.reshape(u_numerical_lexi2, (N,N), order  = "F")