Метод Гаусса с выбором главного элемента

Урок про базовый прямой метод решения СЛАУ — исключение Гаусса — и про то, почему без выбора главного элемента он ненадёжен.

Метод Гаусса приводит систему $Ax = b$ к треугольному виду, последовательно исключая неизвестные (прямой ход), а затем находит их снизу вверх (обратный ход).

Идея: сделать треугольник

Решать треугольную систему легко: в последнем уравнении одна неизвестная — находим её, подставляем в предпоследнее, и так вверх. Поэтому цель прямого хода — обнулить всё под главной диагональю. Берём первое уравнение как «опорное» и вычитаем его из остальных с такими множителями, чтобы в первом столбце под диагональю стали нули. Затем то же со вторым столбцом, используя второе уравнение, и так далее. Получаем верхнетреугольную матрицу — дальше обратный ход.

Зачем выбор главного элемента

Опорный элемент (на который делим, вычисляя множители) называется главным (pivot). Если он окажется нулём — деление невозможно, метод встанет. Если он просто мал — множители огромны, и ошибки округления катастрофически раздуваются. Лекарство — частичный выбор главного элемента: перед каждым шагом меняем строки так, чтобы на диагональ попал наибольший по модулю элемент столбца. Это и спасает от деления на ноль, и резко повышает точность.

def гаусс(A, b):
    n = len(A)
    # расширенная матрица [A | b]
    M = [row[:] + [b[i]] for i, row in enumerate(A)]
    for k in range(n):
        # --- выбор главного элемента: наибольший по модулю в столбце k ---
        p = max(range(k, n), key=lambda i: abs(M[i][k]))
        M[k], M[p] = M[p], M[k]
        # --- исключение ---
        for i in range(k + 1, n):
            f = M[i][k] / M[k][k]
            for j in range(k, n + 1):
                M[i][j] -= f * M[k][j]
    # --- обратный ход ---
    x = [0.0] * n
    for i in range(n - 1, -1, -1):
        s = M[i][n] - sum(M[i][j] * x[j] for j in range(i + 1, n))
        x[i] = s / M[i][i]
    return x

A = [[2.0, 1.0, -1.0],
     [-3.0, -1.0, 2.0],
     [-2.0, 1.0, 2.0]]
b = [8.0, -11.0, -3.0]
x = гаусс(A, b)
print("решение x =", [round(v, 6) for v in x])
# проверка: A*x должно дать b
for i in range(3):
    print(f"уравнение {i+1}: {sum(A[i][j]*x[j] for j in range(3)):.4f}")

Вывод:

решение x = [2.0, 3.0, -1.0]
уравнение 1: 8.0000
уравнение 2: -11.0000
уравнение 3: -3.0000

Сложность и масштаб

Прямой ход стоит примерно $n^3/3$ операций, обратный — $n^2/2$. Итого порядок $O(n^3)$. Для $n = 1000$ это около 300 миллионов операций — доли секунды. Для $n = 1\,000\,000$ (как в инженерных задачах метода конечных элементов) прямой Гаусс уже невозможен — и тут на сцену выходят итерационные методы и учёт разреженности.

Как работает под капотом

Что значит «обнулить элемент» в терминах линейной алгебры? Каждое вычитание строки с множителем — это умножение системы слева на элементарную матрицу. Перестановка строк — умножение на матрицу перестановки $P$. В итоге метод Гаусса с выбором главного элемента эквивалентен разложению $PA = LU$ (об этом — следующий урок): мы неявно строим нижнетреугольную $L$ (из множителей) и верхнетреугольную $U$ (результат прямого хода). Понимание этого превращает Гаусс из «рецепта» в фундамент целого семейства методов.

Частые ошибки

• Гаусс без выбора главного элемента. На «невинной» матрице с малым диагональным элементом получите мусор из-за округления. Частичный pivoting почти бесплатен — используйте всегда.
• Не проверять вырожденность. Если после выбора главный элемент всё равно $\approx 0$, матрица вырождена (нет единственного решения) — нужно это поймать, а не делить на ноль.
• Применять плотный Гаусс к огромной разреженной системе. $O(n^3)$ и заполнение нулей делают это неподъёмным; нужны разреженные или итерационные методы.

Итоги

• Гаусс приводит $Ax=b$ к треугольному виду (прямой ход) и решает его обратным ходом.
• Выбор главного элемента (наибольший по модулю в столбце) спасает от деления на ноль и потери точности.
• Сложность $O(n^3)$ — отлично до тысяч неизвестных, неподъёмно для миллионов.
• По сути это разложение $PA = LU$ — мост к следующему уроку.