The methods for solving systems of linear equations themselves work correctly. Tested with solve. But when trying to solve a large number of matrices, some solution vectors do not converge with the true ones, and therefore a normal distribution for the error does not obtain. The histograms above show that several values ​​are very large - these are incorrectly solved systems of linear equations. Such cases were derived and solved separately by the same functions, but the correct solution was obtained

import numpy as np
import matplotlib.pyplot as plt

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_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)

        b = np.array([1, 1, 1, 1, 1, 1]).astype(np.float64)
     
        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)
  
        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()

I don’t understand where the abnormally large sup_norm and rmse values ​​appear, which prevent a normal distribution from being obtained. As you can see in the histogram, all true values ​​are close to the first column, therefore they accumulate in one place on the graph. I want to eliminate such big mistakes

You didn't correct prevented singular matrix in your random choice.

If you try to plot np.linalg.det(A) vs rmse, you'll see that it is only when determinant is very close to 0 that high values of RMSE occur

Add

determinants.append(np.linalg.det(A))

after the similar lines for errors_rmse and errors_sup_norm,

And then if you plot

plt.scatter(determinants, errors_rmse)

Or in logarithmic scale

Try to restrict the random choosing of matrix to clearly non-singular matrix, and then you'll never see RMSE

The two errors in the way you tried to prevent singular matrix are

Error 1: recursion error

Following code is void:

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

It calls back recursively tridiagonal_matrix, but drops the result. So what you return at the end is just the original A, even if its determinant was 0.

Note that one correction could be

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

But if you think at the recursion nightmare it starts... while and the recursion are a bit redundant here If you want to use recursion you could do it that way

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)
    if np.linalg.det(A) == 0:
        return tridiagonal_matrix(n)
    return A

That would be the elegant way to do it in scheme or caml. But in python, less so. Because python is not a terminal recursion language. Meaning that if it last too long before finding a matrix, it would fail because of stack overflow (the eponymous error)

So it would be better to drop the recursion and keep the while

def tridiagonal_matrix(n):
    while True:
        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)
        if np.linalg.det(A)!=0:
            return A

If you really prefer to avoid the while True: ... if ...: break (some people dogmatically refuse any while True), you can

def tridiagonal_matrix(n):
    A=np.array([[0]])
    while np.linalg.det(A)==0:
        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)
    return A

Error 2 : never compare floats with ==

not being ==0 is not a sufficient condition. Firstly, floating points are not exact numbers. See infamous 0.1+0.2!=0.3.

And secondly, and consequently, even when comparing to exact numbers (like 0), well, because of the accumulation of numeric error, most 0 are never exact 0.

To reuse my previous 0.1+0.2-0.3 example in a linear algebra context np.linalg.det(0.1*np.identity(6)+0.2*np.identity(6)-0.3*np.identity(6)) is not 0.

So you always need a tolerance. And in this case, since you apparently just want to check numerical error of the method (not of the initial matrix) that tolerance has no reason to be very lenient.

So I would

def generate_random_matrix(n = 6 ):
    while True:
        A = (np.random.random((n, n)).astype(np.float64) * 2 - 1).astype(np.float64)
        if not np.isclose(np.linalg.det(A),0, atol=0.1):
            return A

def tridiagonal_matrix(n):
    while True:
        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)
        if not np.isclose(np.linalg.det(A),0, atol=0.1):
            return A

But if that is "too easy" for your need, you can reduce atol

With this way of ensuring matrix are not singular, you get result you expected for "gauss" method

And proof that there is a problem with your thomas implementation (your initial histogram seem to show something that is almost always working. But that is because of the very few cases of huge RMSE due to singular matrices. But if you remove singular matrices, and therefore zoom on the "0" bar, you see that the general case is a always too big RMSE. With the "rmse≈0" being just one random rmse value, no more probable than other)

tl;dr

1. That distribution is due to singular matrices
2. Your way of "retrying" when matrix is singular is ignored, since you always return the first try
3. Your way of testing if a matrix is singluar (with det==0) is too strict. In numerical computation with floats, there is always a tolerance
4. Once singular matrices really avoided, you see that the distribution is centered on very low error value for Gauss. And even the very rare extreme values of this distribution are still very low (order of 10⁻¹⁰ at most)
5. On the other hand, once singular matrices really avoided for Thomas, you see that distribution, sure, doesn't contain crazy values for error, but RMSE>1 is a very typical case, proving simply that your Thomas implemetation is not working (it is not "not working in some strange cases"; it is never working, and the few 0 — and obviously, if you were zooming on only the first bar, you would see that this 0 bar is made of lot of 0.2, 0.3, 0.4, 0.8, ... — the very few "=0 within numerical error margin" are like the 2 times a day a stopped clock is right).

Since this was not your question, and I am a bit lazy, I won't try to understand why your Thomas implementation is not working, but it is not. I guess your initial objective was to see if it were working or not, so mission accomplished :D