Метод Якоби для уравнения Лапласа с граничными условиями Неймана

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[i]) - 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)