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

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

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

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

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

Идея: из $i$-го уравнения выразим $x_i$ через остальные неизвестные: $x_i = (b_i - \sum_{j \ne i} a_{ij} \cdot 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}| \gt \sum_{j \ne i}|a_{ij}|$. Наша матрица ему удовлетворяет ($10 \gt 1+2$ и т.д.), поэтому методы сошлись. Без преобладания итерации могут расходиться — тогда переставляют уравнения, чтобы усилить диагональ, или применяют более мощные методы (сопряжённых градиентов, GMRES).

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

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

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

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

Итоги

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