Add files via upload

poisson_solver_five_point.py

@@ -0,0 +1,120 @@
+import numpy as np
+from scipy.sparse import lil_matrix, csr_matrix
+from scipy.sparse.linalg import spsolve
+import matplotlib.pyplot as plt
+
+# Define the domain and mesh size
+n = 8
+h = 2**(-n)
+x = np.arange(0, 1+h, h)
+y = np.arange(0, 1+h, h)
+X, Y = np.meshgrid(x, y)
+N = len(x)-2  # number of interior nodes
+
+
+# Define the right-hand side function f(x,y) and the boundary conditions g(x,y)
+def f(x, y):
+    return -20*x**4*y**3-12*x**2*y**5 - 17*np.sin(x*y)*(x**2+y**2)
+
+def g(x, y):
+    return x**4*y**5 - 17*np.sin(x*y)
+
+# Create the sparse matrix A using the five-point stencil
+A = lil_matrix((N**2, N**2))
+for i in range(N):
+    for j in range(N):
+        k = i*N + j
+        A[k, k] = -4
+        if i > 0:
+            A[k, k-N] = 1
+        if i < N-1:
+            A[k, k+N] = 1
+        if j > 0:
+            A[k, k-1] = 1
+        if j < N-1:
+            A[k, k+1] = 1
+A = csr_matrix(A)  # convert to compressed sparse row format
+
+# Define the vector b by computing the values of f(x,y) and g(x,y) at each node
+b = np.zeros(N**2)
+for i in range(N):
+    for j in range(N):
+        k = i*N+j
+        b[k] = f(x[i+1], y[j+1])*h**2
+        if i == 0:
+            b[k] -= g(x[i], y[j+1])
+        if i == N-1:
+            b[k] -= g(x[i+2], y[j+1])
+        if j == 0:
+            b[k] -= g(x[i+1], y[j])
+        if j == N-1:
+            b[k] -= g(x[i+1], y[j+2])
+
+# Solve the linear system Ax = b using a sparse direct solver
+u = spsolve(A, b)
+
+# Reshape the solution vector x into a 2D array and plot the solution
+U = np.zeros((N+2, N+2))
+U[1:-1, 1:-1] = u.reshape((N, N))
+
+fig = plt.figure()
+ax = fig.add_subplot(111, projection='3d')
+ax.plot_surface(X, Y, U, cmap='jet')
+ax.set_xlabel('x')
+ax.set_ylabel('y')
+ax.set_zlabel('u')
+ax.set_title('Numerical solution of -∆u=f')
+plt.show()
+
+plt.contourf(X, Y, U, cmap='jet')
+plt.colorbar()
+plt.xlabel('x')
+plt.ylabel('y')
+plt.title('Numerical solution of -∆u=f')
+plt.show()
+
+# Compute exact solution
+G = np.zeros((N+2, N+2))
+for i in range(N+2):
+    for j in range(N+2):
+        if i==0:
+           U[0,j]=g(x[i], y[j])
+        if i==N+1:
+           U[-1,j]=g(x[-1], y[j])
+        if j==0:
+           U[i,0]=g(x[i], y[j])
+        if j==N+1:
+           U[i,-1]=g(x[i], y[j])
+
+# Compute the error norms for two finest meshes
+h1 = 2 ** (-n)  # Mesh size for the first finest mesh
+h2 = 2 ** (-(n+1))  # Mesh size for the second-finest mesh
+
+# Compute the errors for the two finest meshes
+error_inf = np.max(np.abs(U - G))
+error_l2 = np.sqrt(np.sum((U - G)**2))
+
+error_inf2 = np.max(np.abs(U - G))
+error_l22 = np.sqrt(np.sum((U - G)**2))
+
+# Compute the orders of convergence
+order_inf = np.log2(error_inf / error_inf2)
+order_l2 = np.log2(error_l2 / error_l22)
+
+# Print the results
+print("Mesh 1 (h = {})".format(h1))
+print("||u - uh||_inf: ", error_inf)
+print("||u - uh||_l2: ", error_l2)
+print()
+
+print("Mesh 2 (h = {})".format(h2))
+print("||u - uh||_inf: ", error_inf2)
+print("||u - uh||_l2: ", error_l22)
+print()
+
+print("Orders of Convergence:")
+print("Order (||u - uh||_inf): ", order_inf)
+print("Order (||u - uh||_l2): ", order_l2)
+
+
+

poisson_solver_nine_point.py

@@ -0,0 +1,177 @@
+import numpy as np
+from scipy.sparse import lil_matrix, csr_matrix
+from scipy.sparse.linalg import spsolve
+import matplotlib.pyplot as plt
+
+# Define the domain and mesh size
+def nine_point_stencil(n):
+  h = 2**(-n)
+  x = np.arange(0, 1+h, h)
+  y = np.arange(0, 1+h, h)
+  X, Y = np.meshgrid(x, y)
+  N = len(x)-2  # number of interior nodes
+
+
+  # Define the right-hand side function f(x,y) and the boundary conditions g(x,y)
+  def f(x, y):
+      return -20 * x**4 * y**3 - 12 * x**2 * y**5 - 17 *np.sin(x*y) * (x**2 + y**2)
+
+  def g(x, y):
+      return x**4 * y**5 - 17 * np.sin(x*y)
+
+  # Create the sparse matrix A using the five-point stencil
+  gamma = -20                # Main point coefficient
+  alfa = 4                  # Right, Left, Up, and Down neighbors coefficient
+  betta = 1                 # Corner neighbors coefficient
+  A = lil_matrix((N **2, N ** 2))
+  for i in range(N):
+      for j in range(N):
+          k = i * N + j
+          #The main point
+          A[k, k] = gamma
+
+          if j< N - 1:
+              # Right neighbor
+              A[k, k + 1] = alfa
+              # Left neighbor
+              A[k + 1, k ] = alfa
+
+              if i < N - 1:
+                # Up_Left neighbor
+                 A[k + 1, k + N] = betta
+                 # Down_Right neighbor
+                 A[k + N, k + 1] = betta
+
+          if i < N - 1:
+
+              # Up neighbor
+              A[k, k + N] = alfa
+              # Down neighbor
+              A[k + N, k] = alfa
+
+              if j< N - 1:
+                # Up_Right neighbor
+                  A[k, k + N+1] = betta
+                # Down_Left neighbor
+                  A[k + N+1, k] = betta
+
+  A = csr_matrix(A)  # convert to compressed sparse row format
+
+
+  # Define the vector b by computing the values of f(x,y) and g(x,y) at each node
+
+  b = np.zeros(N**2)
+  for i in range(N):
+      for j in range(N):
+          k = i*N+j
+          b[k] = - f(x[i+1], y[j+1])*6*h**2
+
+          if i == 0:
+             if j == 0: 
+                b[k] -= alfa * g(x[0], y[1]) + alfa * g(x[1], y[0]) + betta * g(x[0], y[0]) + betta * g(x[0], y[2]) + betta * g(x[2], y[0])
+
+             if j == N-1: 
+                b[k] -= alfa * g(x[0], y[N]) + alfa * g(x[1], y[N+1]) + betta * g(x[0], y[N+1]) + betta * g(x[0], y[N-1]) + betta * g(x[2], y[N+1])
+
+             if j > 0 and j < N-1:
+                b[k] -= betta * g(x[0], y[j]) + alfa * g(x[0], y[j+1]) + betta * g(x[0], y[j+2])
+
+
+
+          if i == N-1:
+              if j == 0: 
+                b[k] -= alfa * g(x[N+1], y[1]) + alfa * g(x[N], y[0]) + betta * g(x[N-1], y[0]) + betta * g(x[N+1], y[0]) + betta * g(x[N+1], y[2])
+
+              if j == N-1: 
+                b[k] -= alfa * g(x[N+1], y[N]) + alfa * g(x[N], y[N+1]) + betta * g(x[N+1], y[N+1]) + betta * g(x[N-1], y[N+1]) + betta * g(x[N+1], y[N-1])
+
+              if j > 0 and j < N-1:
+                b[k] -= betta * g(x[N+1], y[j]) + alfa * g(x[N+1], y[j+1]) + betta * g(x[N+1], y[j+2])   
+
+
+          if i > 0 and i < N-1:
+             if j == 0:
+                b[k] -= betta * g(x[i], y[0]) + alfa * g(x[i+1], y[0]) + betta * g(x[i+2], y[0])
+
+             if j == N-1:
+                b[k] -= betta * g(x[i], y[N+1]) + alfa * g(x[i+1], y[N+1]) + betta * g(x[i+2], y[N+1])
+
+
+  # Solve the linear system Ax = b using a sparse direct solver
+  u = spsolve(A, b)
+
+  # Reshape the solution vector x into a 2D array and plot the solution
+  U = np.zeros((N+2, N+2))
+  U[1:-1, 1:-1] = u.reshape((N, N))
+
+
+  # Set boundary conditions
+  for i in range(N+2):
+      for j in range(N+2):
+
+          if i == 0:
+             U[0, j]=g(x[i], y[j])
+
+          if i == N+1:
+             U[N+1, j]=g(x[i], y[j])
+
+          if j == 0:
+             U[i, 0]=g(x[i], y[j])
+
+          if j == N+1:
+             U[i, N+1]=g(x[i], y[j])
+
+  u_exact = np.zeros((N+2, N+2))
+  for i in range(N+2):
+      for j in range(N+2):
+        u_exact [i,j]=  g(x[i], y[j]) 
+
+
+  error_inf_matrix = np.abs(U-u_exact)
+  error_inf = np.max(error_inf_matrix)
+  error_l2 = np.sqrt(np.sum((U-u_exact)**2))
+
+
+  fig = plt.figure(1)
+  ax = fig.add_subplot(111, projection='3d')
+  ax.plot_surface(X, Y, U, cmap='jet')
+  ax.set_xlabel('x')
+  ax.set_ylabel('y')
+  ax.set_zlabel('u')
+  ax.set_title('Numerical solution of -∆u=f')
+  plt.show()
+
+  plt.contourf(X, Y, U, cmap='jet')
+  plt.colorbar()
+  plt.xlabel('x')
+  plt.ylabel('y')
+  plt.title('Numerical solution of -∆u=f')
+  plt.show()
+
+  plt.contourf(X, Y, error_inf_matrix, cmap='jet')
+  plt.colorbar()
+  plt.xlabel('x')
+  plt.ylabel('y')
+  plt.title('Error_inf')
+  plt.show()
+
+
+  return error_inf, error_l2
+errors_inf = []
+errors_l2 = []
+i= 0
+for n in [7,8]:
+  i += i
+  result = nine_point_stencil(n)
+  errors_inf.append (result [0])
+  errors_l2 .append (result [1])
+
+odredr_inf = np.log2 (errors_inf[0]/errors_inf[1])
+print (odredr_inf)
+print (errors_inf)
+
+odredr_l2 = np.log2 (errors_l2[0]/errors_l2[1])
+print (odredr_l2)
+print (errors_l2)
+
+