Gauss algorithm

import numpy as np


def gauss_elimination(a, b):
    """
    Решава системата линейни уравнения Ax = b с метод на Гаус.

    Аргументи:
    a -- матрицата на коефициентите (NxN)
    b -- векторът с резултати (Nx1)
    """
    n = len(b)

    # Разширена матрица
    aug_matrix = np.hstack((a, b.reshape(-1, 1)))

    # Привеждане към горнотриъгълна форма
    for i in range(n):
        # Намиране на максималния елемент в колоната за по-голяма стабилност
        max_row = i + np.argmax(abs(aug_matrix[i:, i]))

        # Размяна на текущия ред с този с максимален елемент
        aug_matrix[[i, max_row]] = aug_matrix[[max_row, i]]

        # Нормализация на реда
        aug_matrix[i] = aug_matrix[i] / aug_matrix[i, i]

        # Зануляване на елементите под диагоналния елемент
        for j in range(i + 1, n):
            aug_matrix[j] = aug_matrix[j] - aug_matrix[i] * aug_matrix[j, i]

    # Обратно заместване
    x = np.zeros(n)
    for i in range(n - 1, -1, -1):
        x[i] = aug_matrix[i, -1] - np.sum(aug_matrix[i, i + 1:n] * x[i + 1:n])

    return x


# Пример
a = np.array([
    [1, 1, 1],
    [2, 3, 5],
    [4, 0, 5]], dtype=float)

b = np.array([5, 8, 2], dtype=float)

solution = gauss_elimination(a, b)
print("Решението x е:", solution)




def gauss_elimination(A, b):
    """
    Решава система линейни уравнения Ax = b използвайки метода на Гаус

    Параметри:
    A (numpy.array): Матрица на коефициентите
    b (numpy.array): Вектор на свободните членове

    Връща:
    numpy.array: Решение на системата
    """
    # Превръщаме входните данни в numpy масиви
    A = np.array(A, dtype=np.float64)
    b = np.array(b, dtype=np.float64)
    n = len(A)

    # Разширена матрица [A|b]
    Ab = np.concatenate((A, b.reshape(n, 1)), axis=1)

    # Права стъпка (елиминация)
    for i in range(n):
        # Намираме максималния елемент в текущата колона за частично пивотиране
        pivot_row = i + np.argmax(abs(Ab[i:, i]))
        if pivot_row != i:
            Ab[i], Ab[pivot_row] = Ab[pivot_row].copy(), Ab[i].copy()

        pivot = Ab[i, i]
        if abs(pivot) < 1e-10:
            raise ValueError("Матрицата е сингулярна или почти сингулярна")

        # Нормализиране на текущия ред
        Ab[i] = Ab[i] / pivot

        # Елиминиране на елементите под главния диагонал
        for j in range(i + 1, n):
            multiplier = Ab[j, i]
            Ab[j] = Ab[j] - multiplier * Ab[i]

    # Обратна стъпка (обратна субституция)
    x = np.zeros(n)
    for i in range(n - 1, -1, -1):
        x[i] = Ab[i, -1] - np.sum(Ab[i, i + 1:n] * x[i + 1:])

    return x


# Пример за използване:
def example_usage():
    # Система:
    # 2x + y - z = 8
    # -3x - y + 2z = -11
    # -2x + y + 2z = -3

    A = [
        [1, 1, 1],
        [2, 3, 5],
        [4, 0, 5]
    ]
    b = [5, 8, 2]

    try:
        solution = gauss_elimination(A, b)
        print("Решение:", solution)

        # Проверка на решението
        print("\nПроверка (Ax = b):")
        print("Ax =", np.dot(A, solution))
        print("b =", b)

    except ValueError as e:
        print("Грешка:", e)


if __name__ == "__main__":
    example_usage()