From 5baf42bcdf8f722c5435060aabee90be89197914 Mon Sep 17 00:00:00 2001
From: javak87 <javak87@mi.fu-berlin.de>
Date: Tue, 4 May 2021 16:16:58 +0430
Subject: [PATCH] first edition FDM
---
Poisson_equation_FDM.py | 88 ++++++++++++++++++++---------------------
1 file changed, 44 insertions(+), 44 deletions(-)
diff --git a/Poisson_equation_FDM.py b/Poisson_equation_FDM.py
index 4956c8f..b9f5470 100644
--- a/Poisson_equation_FDM.py
+++ b/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()
+
--
GitLab