[Python] Jacobi method converges faster than Gauss-Seidel in case of Laplace equation

It is need to solve the problem using the Jacobi method: d^T/dx^2 + d^2T/dy^2 = 0 with conditions at the boundaries of the square: T(0, y) = 1, T(x, 1) = 1, dT/dx(1, y) = 1 - y, dT/dy(x, 0) = x.

Conditions for stopping iterations: max|T_new - T| < d_t and max|R_{i, j}| < d_R, where R is the residual matrix.

I approximate the boundary conditions on derivatives by 2-order finite-difference schemes. By the condition of the problem, the step along the x and y axes is h = 1/20. My implementation of the numerical method results in 645 iterations for convergence in the case of the Jacobi method and 999 in the case of Gauss-Seidel, which is most likely incorrect. I think that the problem is in the Jacobi method. Please help me find out why this can happen. I enclose my code for the Jacobi method.

h = 1/20
grid = np.arange(0, 1+h, h)
N = len(grid)

def jacobi_method(M):
T = M.copy()
T[0, :] = 1
T[:, -1] = 1
T_prev = T.copy()
for i in range(1, N - 1):
for j in range(1, N - 1):
T[i, j] = ( T_prev[i+1, j] + T_prev[i-1, j] + T_prev[i, j+1] + T_prev[i, j-1]) /4

T[-1, :] = (4*T[-2, : ] - T[-3, :] + 2*h*(1- grid))/3
T[:,0] = (4*T[:, 1] - T[:, 2] - 2*h*grid)/3

return T

def residual(T):
residual = np.zeros((len(T)-2, len(T)-2))
for i in range(1, len(T) - 1):
for j in range(1, len(T) - 1):
residual[i-1, j-1] = ((T[i+1, j] - 2*T[i, j] + T[i-1, j]) / (h**2) +
(T[i, j+1] - 2*T[i, j] + T[i, j-1]) / (h**2))

return residual

T1 = np.ones((N, N))
dT = 1
deltaT = 10e-7
R = 1
deltaR = 10e-5
iterations = []
dTs = []
Rs = []
i = 0
while dT >= deltaT or R >= deltaR :
T_n = jacobi_method(T1)
dT = np.max(abs(T_n - T1))
dTs.append(dT)
R = np.max(abs(residual(T_n)))
Rs.append(R)
T1 = T_n.copy()
i += 1
iterations.append(i)

Continue reading...