Итерационные методы: Якоби и Зейдель

Урок про итерационное решение СЛАУ методами Якоби и Зейделя: когда оно лучше прямого и при каком условии сходится.

Метод Якоби выражает каждое неизвестное из своего уравнения через остальные и пересчитывает все одновременно. Метод Зейделя — то же, но сразу использует уже обновлённые значения, ускоряя сходимость.

Зачем итерации, если есть Гаусс

Прямой Гаусс точен, но стоит O(n³) и для n = 10^6 неподъёмен. А такие системы реальны: дискретизация уравнений в физике даёт миллионы неизвестных, но матрица разреженная — в каждой строке лишь несколько ненулевых элементов. Итерационные методы не трогают нули: один шаг стоит O(число ненулей), и при хорошей сходимости несколько десятков шагов дают нужную точность дешевле, чем кубический Гаусс. Плюс память: храним только ненулевые элементы.

Идея: из i-го уравнения выразим x_i через остальные неизвестные: x_i = (b_i − Σ_{j≠i} a_{ij}·x_j) / a_{ii}. Это формула пересчёта. Стартуем с нулей и крутим.

def якоби(A, b, итераций=15):
    n = len(A)
    x = [0.0] * n
    for it in range(1, итераций + 1):
        x_new = [0.0] * n
        for i in range(n):
            s = sum(A[i][j] * x[j] for j in range(n) if j != i)
            x_new[i] = (b[i] - s) / A[i][i]
        x = x_new            # все компоненты обновляются разом
        if it % 3 == 0:
            print(f"Якоби, итерация {it:2}: {[round(v, 5) for v in x]}")
    return x

A = [[10.0, -1.0, 2.0],
     [-1.0, 11.0, -1.0],
     [2.0, -1.0, 10.0]]
b = [6.0, 25.0, -11.0]
якоби(A, b)

Вывод:

Якоби, итерация  3: [1.02127, 2.27769, -1.08673]
Якоби, итерация  6: [1.04346, 2.26906, -1.08141]
Якоби, итерация  9: [1.04326, 2.26923, -1.08174]
Якоби, итерация 12: [1.04327, 2.26923, -1.08173]
Якоби, итерация 15: [1.04327, 2.26923, -1.08173]

Зейдель: используем свежие значения

Якоби обновляет все x_i одновременно по «старым» значениям. Но при пересчёте x_2 мы уже знаем новое x_1 — почему бы не использовать его сразу? Это и есть метод Зейделя: обновлённые компоненты немедленно идут в дело. Обычно он сходится примерно вдвое быстрее Якоби.

def зейдель(A, b, итераций=10):
    n = len(A)
    x = [0.0] * n
    for it in range(1, итераций + 1):
        for i in range(n):
            s = sum(A[i][j] * x[j] for j in range(n) if j != i)
            x[i] = (b[i] - s) / A[i][i]   # сразу пишем в x — свежие значения
        if it % 2 == 0:
            print(f"Зейдель, итерация {it:2}: {[round(v, 6) for v in x]}")
    return x

A = [[10.0, -1.0, 2.0], [-1.0, 11.0, -1.0], [2.0, -1.0, 10.0]]
b = [6.0, 25.0, -11.0]
зейдель(A, b)

Вывод:

Зейдель, итерация  2: [1.030182, 2.276628, -1.078374]
Зейдель, итерация  4: [1.043297, 2.269235, -1.081736]
Зейдель, итерация  6: [1.043269, 2.269231, -1.081731]
Зейдель, итерация  8: [1.043269, 2.269231, -1.081731]
Зейдель, итерация 10: [1.043269, 2.269231, -1.081731]

Зейдель «устаканился» к 6-й итерации, Якоби — к 12-й: вдвое быстрее, как и обещано.

Условие сходимости

Оба метода сходятся не всегда. Достаточное условие — диагональное преобладание: в каждой строке модуль диагонального элемента больше суммы модулей остальных, |a_{ii}| > Σ_{j≠i}|a_{ij}|. Наша матрица ему удовлетворяет (10 > 1+2 и т.д.), поэтому методы сошлись. Без преобладания итерации могут расходиться — тогда переставляют уравнения, чтобы усилить диагональ, или применяют более мощные методы (сопряжённых градиентов, GMRES).

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

Любой такой метод записывается как x_{n+1} = T·x_n + c, где T — матрица перехода. Сходимость определяется её спектральным радиусом (наибольшим по модулю собственным значением): процесс сходится тогда и только тогда, когда ρ(T) < 1 — прямой аналог условия |g'| < 1 из простой итерации. Чем меньше ρ(T), тем быстрее сходимость. Диагональное преобладание — простое достаточное условие, гарантирующее ρ(T) < 1.

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

Запускать без диагонального преобладания. Можно получить расходимость; сначала переставьте строки/столбцы ради сильной диагонали.
Останавливаться по числу итераций, а не по невязке. Контролируйте ||Ax − b|| или ||x_{n+1} − x_n|| — иначе либо переплатите, либо недосчитаете.
Применять плотную реализацию к разреженной системе. Весь выигрыш — в пропуске нулей; используйте разреженное хранение.

Итоги

Итерационные методы уточняют приближение, не трогая нули — это путь к большим разреженным системам.
Якоби обновляет все компоненты разом; Зейдель сразу использует свежие — и сходится примерно вдвое быстрее.
Достаточное условие сходимости — диагональное преобладание (формально ρ(T) < 1).
Останавливаться следует по невязке/изменению, а не по фиксированному числу шагов.