Метод Ричардсона

У меня такая проблема, если я прошу сделать print(Richardson(2)), то всё работает, но если я поменяю 2 на 12, то мне выдаёт в этом методе ошибку ERROR thonny.backend: PROBLEM WITH THONNY'S BACK-END. И если я поставлю w = 1.99/h_norm (как на самом деле должно быть), то ни в одном случае метод не работает. Можете подсказать решение, как реализовать метод Ричардсона для решения системы уравнений? Есть вариант, что этот метод нужно сделать циклом, но я не понимаю как.

import numpy as np

def h(n):
    maatriks = []
    for i in range(1, n+1):
        hi = []
        for j in range(1, n+1):
            h_i_j = 1/(i+j-1)
            hi.append(h_i_j)
        maatriks.append(hi)
    return maatriks

def x(n):
    xi = []
    for i in range(0, n):
        xi.append(1)
    return xi

def b(n):
    bi = []
    h1 = h(n)
    x1 = x(n)
    b_i = np.dot(h1, x1)
    for i in b_i:
        bi.append(i)
    return bi

def b_theta(n):
    bi_theta = []
    h1 = h(n)
    x1 = x(n)
    b_i_theta = np.dot(h1, x1)
    for i in b_i_theta:
        bi_theta.append(round(i, 9))
    return bi_theta

def Gauss(n):
    h1 = h(n)
    b1 = b(n)
    b_theta1 = b_theta(n)
    x1_1 = np.linalg.solve(h1, b1)
    x2_1 = np.linalg.solve(h1, b_theta1)
    x_tärn = x(n)
    x1_1_1 = x1_1 - x_tärn
    x2_1_1 = x2_1 - x_tärn
    x1_norm = np.linalg.norm(x1_1_1, ord = np.inf)
    x2_norm = np.linalg.norm(x2_1_1, ord = np.inf)
    return x1_norm, x2_norm

def pöörd(n):
    h1 = h(n)
    b1 = b(n)
    b_theta1 = b_theta(n)
    h_pöörd = np.linalg.inv(h1)
    x1_1 = np.dot(h_pöörd, b1)
    x2_1 = np.dot(h_pöörd, b_theta1)
    x_tärn = x(n)
    x1_1_1 = x1_1 - x_tärn
    x2_1_1 = x2_1 - x_tärn
    x1_norm = np.linalg.norm(x1_1_1, ord = np.inf)
    x2_norm = np.linalg.norm(x2_1_1, ord = np.inf)
    return x1_norm, x2_norm

def D(H):
    Di = []
    for i in H:
        d = []
        for j in i:
            if H.index(i) == i.index(j):
               d.append(j)
            else:
                d.append(0)
        Di.append(d)
    return Di

def R(H, D):
    M = []
    for i, j in zip(H, D):
        m = []
        for k, l in zip(i, j):
            s = k - l
            m.append(s)
        M. append(m)
    return M

def jacobi(x0, b, maxit):
    if maxit == 0:
        return x0
    else:
        H = h(len(x0))
        D1 = D(H)
        R1 = R(H, D1)
        D_poord = np.linalg.inv(D1)
        x_1 = np.dot(D_poord, b) - np.dot(np.dot(D_poord, R1), x0)
        maxit -= 1
        return jacobi(x_1, b, maxit)

def xm(w, x0, x0_theta, H, b, b_theta, theta, n):
    xm_1 = x0 - w*(np.dot(H, x0) - b)
    xm_1_theta = x0 - w*(np.dot(H, x0_theta) - b_theta)
    norm1 = np.dot(H, xm_1) - b
    norm2 = np.dot(H, xm_1_theta) - b_theta
    if np.linalg.norm(norm1, ord = 2) <= theta*np.sqrt(n) or np.linalg.norm(norm2, ord = 2) <= theta*np.sqrt(n):
        return xm_1, xm_1_theta
    else:
        return xm(w, xm_1, xm_1_theta, H, b, b_theta, theta, n)

def Richardson(n):
    h1 = h(n)
    b1 = b(n)
    b_theta1 = b_theta(n)
    h_norm = np.linalg.norm(h1, ord = 2)
    w = 1
    x1 = xm(w, 0, 0, h1, b1, b_theta1, 1e-9, n)
    x1_1 = x1[0]
    x2_1 = x1[1]
    x_tärn = x(n)
    x1_1_1 = x1_1 - x_tärn
    x2_1_1 = x2_1 - x_tärn
    x1_norm = np.linalg.norm(x1_1_1, ord = np.inf)
    x2_norm = np.linalg.norm(x2_1_1, ord = np.inf)
    return x1_norm, x2_norm

print(Gauss(12))
print(pöörd(12))
print(Richardson(2))