first edition FDM

Poisson_equation_FDM.py

@@ -17,6 +17,7 @@ class FivePointStencil_2D:
     '''
 
     def __init__ (grid_number_x:'int',grid_number_y:'int', boundary_x:'tuple', boundary_y:'tuple'):
+        
         '''
         Parameters:
                     grid_number_x (int): the number of grid in the x direction
@@ -30,50 +31,49 @@ class FivePointStencil_2D:
         self.boundary_x = boundary_x
         self.boundary_y = boundary_y
 
-        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)
-
-            #xx, yy = np.meshgrid(grid, grid, sparse=True)
-
-            '''
-            initialize boundary condition
-            '''
-
-            '''
-            # 100*100 the problem dimension
-
-            # construct A matrxi (100^2*100^2)
-            eins = np.ones((3,))
-            A = -2*np.eye(3)+np.diag(eins, k=1)[0:-1, 0:-1]+np.diag(eins, k=-1)[0:-1, 0:-1]
-
-            new = np.hstack((base,np.eye(3) ))
-            new =np.hstack((new,np.eye(3) ))
-            new2= np.hstack((np.eye(3), base ))
-            new2= np.hstack((new2, np.eye(3) ))
-            new3 = np.hstack((np.eye(3), np.eye(3) ))
-            new3 = np.hstack((new3, base ))
-            A = np.vstack((new,new2))
-            A = np.vstack((A,new3))
-            '''
-            '''
-            eins = np.ones((9,))
-            A1 = 4*np.eye(9)
-            A11=np.diag(eins, k=1)[0:-1, 0:-1]
-            A12 = np.diag(eins, k=-1)[0:-1, 0:-1]
-            A2 = np.diag(eins, k=3)[0:-3, 0:-3]
-            A3 = np.diag(eins, k=-3)[0:-3, 0:-3]
-            # ADD perodic
-            Aper1= np.diag(eins, k=6)[0:-6, 0:-6]
-            Aper2= np.diag(eins, k=-6)[0:-6, 0:-6]
-            A = A1+A11+A12+A2+A3+Aper1+Aper2
-
-            b= np.array([0.0625, 0.0625, 0.0625])
-            x = np.linalg.solve(A, b)
-            '''
+    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)
+
+        xx, yy = np.meshgrid(x_grid, y_grid, sparse=True)
+        
+        # Initialize all nodes (inner and outer) by the given function
+        values = (np.power(xx, 4))*(np.power(yy, 5)) + np.sin(xx*yy)
+
+        # initialize all inner node by zero
+        values  = values[1:-1, 1:-1]
+        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
+
+        # 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))
+
+        half_matrix_upper = np.hstack((block_matrix,corner_matrix))
+        half_matrix_lower = np.hstack((corner_matrix,block_matrix))
+        A = np.vstack((half_matrix_upper,half_matrix_lower))
+
+        u = np.linalg.solve(A, values)
+
+        return u
+
+if __name__=="__main__":
+
+    fdm_obj = FivePointStencil_2D(grid_number_x=10,grid_number_y=10, boundary_x=[0,1], boundary_y=[0,1])
+    solution = fdm_obj.boundary_initialization()
+