scipy.linalg: решение СЛАУ методом Гаусса

Система линейных уравнений — фундамент инженерных расчётов. Решим её «руками» методом Гаусса и одной строкой в SciPy.

СЛАУ (система линейных алгебраических уравнений) записывается как A·x = b, где A — матрица коэффициентов, b — вектор правых частей, x — искомые неизвестные.

Зачем это нужно

К решению СЛАУ сводится огромное число задач: баланс химических реакций, расчёт электрических цепей (законы Кирхгофа), метод наименьших квадратов, конечно-элементный анализ конструкций, экономические модели «затраты-выпуск». Если вы умеете решать A·x = b — вы держите в руках половину прикладной математики.

Пример системы

Пусть надо решить:

2x + y  = 5
x  + 3y = 10

В матричной форме: A = [[2, 1], [1, 3]], b = [5, 10]. Ответ — x = 1, y = 3 (проверьте подстановкой: 2·1 + 3 = 5, 1 + 3·3 = 10).

Идея метода Гаусса

Метод Гаусса — это аккуратное «обнуление» переменных. Прямой ход: вычитаем уравнения друг из друга, чтобы под главной диагональю остались нули (получаем «ступенчатую» матрицу). Обратный ход: из последнего уравнения находим последнюю переменную, подставляем выше — и так до первой.

Реализация на чистом Python

Вот метод Гаусса с обратной подстановкой на stdlib — он реально решает нашу систему:

def solve(A, b):
    n = len(A)
    # копии, чтобы не портить вход
    M = [row[:] for row in A]
    v = b[:]
    # --- прямой ход: зануляем под диагональю ---
    for col in range(n):
        pivot = M[col][col]
        for row in range(col + 1, n):
            factor = M[row][col] / pivot
            for k in range(col, n):
                M[row][k] -= factor * M[col][k]
            v[row] -= factor * v[col]
    # --- обратный ход ---
    x = [0.0] * n
    for i in range(n - 1, -1, -1):
        s = v[i]
        for j in range(i + 1, n):
            s -= M[i][j] * x[j]
        x[i] = s / M[i][i]
    return x

A = [[2.0, 1.0],
     [1.0, 3.0]]
b = [5.0, 10.0]
print("x, y =", solve(A, b))

Вывод:

x, y = [1.0, 3.0]

А в SciPy — одна строка

То, что мы расписали на двадцать строк, SciPy делает так (код не исполняется в браузере, помечен текстом):

import numpy as np
from scipy.linalg import solve

A = np.array([[2, 1], [1, 3]], dtype=float)
b = np.array([5, 10], dtype=float)
x = solve(A, b)
print(x)   # [1. 3.]

scipy.linalg.solve внутри использует LU-разложение с выбором главного элемента (об этом — следующий урок) — куда устойчивее наивного Гаусса.

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

Наш наивный код делит на M[col][col] — главный элемент (pivot). Если он окажется нулём или очень маленьким, мы получим деление на ноль или катастрофическую потерю точности. Поэтому в настоящих библиотеках (LAPACK, который зовёт SciPy) применяют частичный выбор главного элемента: перед каждым шагом строки переставляют так, чтобы на диагонали оказался максимальный по модулю элемент. Это и есть разница между учебным алгоритмом и промышленным: правильность в краевых случаях.

Сложность метода Гаусса — O(n³). Для системы 1000×1000 это миллиард операций; вот почему делать их в скомпилированном LAPACK, а не в Python-цикле, критично.

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

Нулевой pivot. Наивный Гаусс падает, если на диагонали ноль; нужен выбор главного элемента.
Вырожденная матрица. Если уравнения линейно зависимы (det = 0), решения нет или их бесконечно много — SciPy выдаст ошибку или предупреждение.
Целочисленные массивы. Если не указать dtype=float, деление может округлиться и испортить ответ.

Итог

СЛАУ A·x = b — основа множества инженерных задач.
Метод Гаусса: прямой ход (зануление под диагональю) + обратная подстановка.
«Руками» — для понимания; в проде — scipy.linalg.solve (LAPACK, выбор главного элемента).
Сложность O(n³) — большие системы решают только в скомпилированном коде.