Initial commit

Poisson_equation_FDM.py

+# -*- coding: utf-8 -*-
+
+'''
+solve possion eq in two D
+The right-hand side and the boundary conditions should be chosen such that:
+    
+u(x; y) = (x^4)*(y^5) + 17*sin(xy)
+  
+'''
+import numpy as np
+import matplotlib.pyplot as plt
+
+class FivePointStencil_2D:
+    '''
+    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'):
+        '''
+        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]
+        '''
+        
+        self.grid_number_x = grid_number_x
+        self.grid_number_y = grid_number_y
+        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)
+            '''
+
+
+