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

我不明白异常大的sup_NORAME和RMSE值​​出现在哪里，它们阻止了正态分布的获得.正如您在直方图中看到的那样，所有真值​​都接近第一列，因此它们累积在图形上的一个位置.我想要消除这样的大错误