From 3e14aece473b95937c77afef2dee2fd141c62d3a Mon Sep 17 00:00:00 2001 From: Saeid Zarifikoliaee Date: Fri, 2 Jun 2023 21:32:27 +0200 Subject: [PATCH] Add files via upload --- poisson_solver_five_point.py | 120 ++++++++++++++++++++++++ poisson_solver_nine_point.py | 177 +++++++++++++++++++++++++++++++++++ 2 files changed, 297 insertions(+) create mode 100644 poisson_solver_five_point.py create mode 100644 poisson_solver_nine_point.py diff --git a/poisson_solver_five_point.py b/poisson_solver_five_point.py new file mode 100644 index 0000000..66bb7ec --- /dev/null +++ b/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) + + + diff --git a/poisson_solver_nine_point.py b/poisson_solver_nine_point.py new file mode 100644 index 0000000..73be68a --- /dev/null +++ b/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) + +