second revision still matrix mismach

Poisson_equation_FDM.py

@@ -10,34 +10,31 @@ u(x; y) = (x^4)*(y^5) + 17*sin(xy)
 import numpy as np
 import matplotlib.pyplot as plt
 
-class FivePointStencil_2D:
+class FivePointStencil:
     '''
     This class solve the Poisson equation numerically in 2D
     using the five-point stencil finite difference method.
     '''
 
-    def __init__ (grid_number_x:'int',grid_number_y:'int', boundary_x:'tuple', boundary_y:'tuple'):
+    def __init__ (self,nx,ny):
         
         '''
         Parameters:
-                    grid_number_x (int): the number of grid in the x direction
-                    grid_number_y (int): the number of grid in the y direction
-                    boundary_x (tuple): the boundary condition in x direction like [lower_bound,upper_bound]
-                    boundary_y (tuple): the boundary condition in y direction like [lower_bound,upper_bound]
+                    nx (int): the number of grid in the x direction
+                    ny (int): the number of grid in the y direction
+
         '''
         
-        self.grid_number_x = grid_number_x
-        self.grid_number_y = grid_number_y
-        self.boundary_x = boundary_x
-        self.boundary_y = boundary_y
+        self.nx = nx
+        self.ny = ny
 
     def boundary_initialization (self):
         '''
         Initialization of boundary condition based on the right hand side function
         right hand side function: u(x; y) = (x^4)*(y^5) + 17*sin(xy)
         '''
-        x_grid = np.linspace(self.boundary_x[0], self.boundary_x[1], self.grid_number_x+1)
-        y_grid = np.linspace(self.boundary_y[0], self.boundary_y[1], self.grid_number_y+1)
+        x_grid = np.linspace(0, 1, self.nx+1)
+        y_grid = np.linspace(0, 1, self.ny+1)
 
         xx, yy = np.meshgrid(x_grid, y_grid, sparse=True)
         
@@ -49,17 +46,21 @@ class FivePointStencil_2D:
         values = values.flatten()
 
         # compute hx and hy
-        hx = (self.boundary_x[1] - self.boundary_x[0])/self.grid_number_x
-        hy = (self.boundary_y[1] - self.boundary_y[0])/self.grid_number_y
+        hx =1/self.nx
+        hy = 1/self.ny
 
         # Au - Bg = f
         # compute "A" Matrix   
                 
         diag_coeff = 2*((1/hx**2) + (1/hy**2))
-        eins = (1/hx**2)*np.ones((self.grid_number_x,))
-        block_matrix = diag_coeff * np.eye((self.grid_number_x-1)*(self.grid_number_x-1)) + np.diag(-1*eins, k=1)[0:-1, 0:-1] + np.diag(eins, k=-1)[0:-1, 0:-1]
-
-        corner_matrix = (1/hy**2)*np.eye((self.grid_number_x-1)*(self.grid_number_x-1))
+        eins = (1/hx**2)*np.ones((self.nx+1,))
+        block_matrix = diag_coeff * np.eye((self.nx-1)) 
+        off_diag_upper = np.diag(-1*eins, k=1)
+        off_diag_lower = np.diag(eins, k=-1)
+        
+        block_matrix = block_matrix + off_diag_upper + off_diag_lower
+        
+        corner_matrix = (1/hy**2)*np.eye((self.nx-1)*(self.ny-1))
 
         half_matrix_upper = np.hstack((block_matrix,corner_matrix))
         half_matrix_lower = np.hstack((corner_matrix,block_matrix))
@@ -71,7 +72,9 @@ class FivePointStencil_2D:
 
 if __name__=="__main__":
 
-    fdm_obj = FivePointStencil_2D(grid_number_x=10,grid_number_y=10, boundary_x=[0,1], boundary_y=[0,1])
+    nx=4
+    ny=3
+    fdm_obj = FivePointStencil(nx,ny)
     solution = fdm_obj.boundary_initialization()