Не получается нормальное распределение для вычислительной ошибки методов решения СЛАУ Гаусом и прогонкой

Сами методы решения СЛАУ работают корректно. Проверено при помощи slove. Но при попытке решить большое количество матриц некоторые вектора решений не сходятся с истинными, всвязи с чем не получается нормального распределения для погрешности. На приведенных гистограммах видно, что несколько значений очень большие - это и есть неверно решенные СЛАУ. Такие случаи были выведены и решены отдельно этими же функциями, однако решение было получено верное

import matplotlib.pyplot as plt
from datetime import datetime
from numpy.linalg import norm, solve
import typing as t


def calculate_relative_error(true_solution, computed_solution):
    n = len(true_solution)

    rmse = np.sqrt(np.sum((true_solution - computed_solution)**2))/n
    sup_norm = np.max(np.abs(true_solution - computed_solution))
    return rmse, sup_norm


def generate_symmetric_matrix(n):

    A = (np.random.random((n, n)).astype(np.float64) * 2 - 1).astype(np.float64)

    A = np.dot(A, A.T)
    A += np.eye(n, dtype=np.float64)

    while np.linalg.det(A) == 0:
        generate_symmetric_matrix(n)

    return A

def generate_random_matrix(n = 6 ):

    A = (np.random.random((n, n)).astype(np.float64) * 2 - 1).astype(np.float64)

    while np.linalg.det(A) == 0:
        generate_random_matrix(n)

    return A

def tridiagonal_matrix(n):

    main_diag = (np.random.random(n).astype(np.float64) * 2 - 1).astype(np.float64)
    sub_diag = (np.random.random(n-1).astype(np.float64) * 2 - 1).astype(np.float64)
    super_diag = (np.random.random(n-1).astype(np.float64) * 2 - 1).astype(np.float64)

    A = np.diag(main_diag) + np.diag(sub_diag, k=-1) + np.diag(super_diag, k=1)

    while np.linalg.det(A) == 0:
        tridiagonal_matrix(n=6)

    return A

def gauss(A, b, pivoting):
    n = len(b)

    a = np.hstack((A, b[:, np.newaxis])).astype(np.float64)

    for i in range(n):
        if pivoting:
            max_index = np.argmax(np.abs(a[i:, i])) + i
            a[[i, max_index]] = a[[max_index, i]]

        for j in range(i + 1, n):
            factor = a[j, i] / a[i, i]
            a[j, i:] -= factor * a[i, i:]

    x = np.zeros(n, dtype=np.float64)
    for i in range(n - 1, -1, -1):
        x[i] = (a[i, -1] - np.dot(a[i, i + 1:-1], x[i + 1:])) / a[i, i]
    return x


def thomas(A, b):
    gamma = [-A[0][1] / A[0][0]]
    beta = [b[0] / A[0][0]]

    n = len(b)
    x = [0] * n

    for i in range(1, n):
        if i != n - 1:
            gamma.append(-A[i][i + 1] / (A[i][i - 1] * gamma[i - 1] + A[i][i]))
        beta.append((b[i] - A[i][i - 1] * beta[i - 1]) / (A[i][i - 1] * gamma[i - 1] + A[i][i]))

    x[n - 1] = beta[n - 1]

    for i in range(n - 2, -1, -1):
        x[i] = gamma[i] * x[i + 1] + beta[i]

    return x






num_matrices = 1000
methods = ["gauss_no_pivoting", "thomas"] 


for method in methods:

    errors_rmse = []
    errors_sup_norm = []

    for _ in range(num_matrices):
        A_random = generate_random_matrix(6)
        A = A_random.copy()
        A_tridiagonal = tridiagonal_matrix(6)
        A_symmetric = generate_symmetric_matrix(6)

        b = np.array([1, 1, 1, 1, 1, 1]).astype(np.float64)
        x = np.linalg.solve(A_random, b)
        true_solution = gauss(A_random, b, pivoting=True)

        computed_solution = None
        if method == "gauss_no_pivoting":
            computed_solution = gauss(A_random.copy(), b.copy(), pivoting=False)
        elif method == "thomas":
            computed_solution = thomas(A_tridiagonal.copy(), b.copy())
       
        rmse, sup_norm = calculate_relative_error(true_solution, computed_solution)
        # if rmse <= 1e-3:
        #     print(A)
        #     print(true_solution)
        #     print(computed_solution)
        #     print(rmse)
        #     print(sup_norm)
        errors_rmse.append(rmse)
        errors_sup_norm.append(sup_norm)




    plt.figure(figsize=(10, 5))
    plt.subplot(1, 2, 1)
    plt.hist(errors_rmse, bins=20, color='blue', edgecolor='black')
    plt.title(f'{method} - RMSE Histogram')
    plt.xlabel('RMSE')
    plt.ylabel('Frequency')

    plt.subplot(1, 2, 2)
    plt.hist(errors_sup_norm, bins=20, color='green', edgecolor='black')
    plt.title(f'{method} - Sup Norm Histogram')
    plt.xlabel('Sup Norm')
    plt.ylabel('Frequency')

    plt.tight_layout()
    plt.show()