Метод Якоби решение СЛАУ: формулы, сходимость, пример

Метод Якоби - простейший итерационный способ решения системы линейных алгебраических уравнений $Ax = b$. Вместо прямого исключения переменных он строит последовательность приближений $x^{(0)}, x^{(1)}, x^{(2)}, \ldots$, сходящуюся к точному решению. Метод незаменим, когда матрица $A$ большая и разреженная: дискретизация уравнений теплопроводности, метод конечных разностей, задачи электростатики - там прямое LU-разложение слишком дорого по памяти, а каждая итерация Якоби стоит лишь одно умножение матрицы на вектор.

Постановка задачи

Дана система $Ax = b$ из $n$ уравнений с $n$ неизвестными, где $A \in \mathbb{R}^{n \times n}$, $b \in \mathbb{R}^n$. Предполагаем, что все диагональные элементы ненулевые: $a_{ii} \neq 0$ для каждого $i = 1, \ldots, n$. Если это не так, уравнения переставляют, чтобы на диагонали оказались ненулевые (а лучше - наибольшие по модулю) коэффициенты. Итерация стартует с произвольного начального вектора $x^{(0)}$, чаще всего нулевого, и на каждом шаге пересчитывает все компоненты сразу.

Идея метода Якоби

Выразим из $i$-го уравнения переменную $x_i$ через все остальные:

$$x_i = \frac{1}{a_{ii}}\left(b_i - \sum_{j \neq i} a_{ij} x_j\right).$$

Если подставить в правую часть значения с предыдущего шага, получим итерационную формулу метода Якоби:

$$x_i^{(k+1)} = \frac{1}{a_{ii}}\left(b_i - \sum_{j \neq i} a_{ij} x_j^{(k)}\right), \quad i = 1, 2, \ldots, n.$$

Принципиальная особенность: все компоненты нового вектора $x^{(k+1)}$ вычисляются на основе одного и того же «старого» вектора $x^{(k)}$. Порядок обхода $i$ неважен, и формула легко распараллеливается - каждая компонента считается независимо. Это и отличает Якоби от метода Гаусса–Зейделя, где уже обновлённые на текущем шаге значения сразу идут в дело.

Матричная форма и матрица перехода

Расщепим $A$ на сумму трёх матриц: $A = D + L + U$, где $D$ - диагональ, $L$ - строго нижнетреугольная часть, $U$ - строго верхнетреугольная. Тогда одна итерация Якоби записывается как $D x^{(k+1)} = b - (L + U)x^{(k)}$, откуда

$$x^{(k+1)} = -D^{-1}(L + U)\, x^{(k)} + D^{-1} b = T_J\, x^{(k)} + c.$$

Матрица $T_J = -D^{-1}(L + U)$ называется матрицей перехода (итерационной матрицей) метода Якоби, а вектор $c = D^{-1} b$. Это запись стандартного вида $x^{(k+1)} = T x^{(k)} + c$, к которому сводятся все простые итерационные методы. Поведение метода целиком определяется свойствами матрицы $T_J$: точное решение $x^{*}$ является её неподвижной точкой ($x^{*} = T_J x^{*} + c$), а ошибка $e^{(k)} = x^{(k)} - x^{*}$ эволюционирует как $e^{(k+1)} = T_J\, e^{(k)}$, то есть $e^{(k)} = T_J^{\,k}\, e^{(0)}$.

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

Из соотношения $e^{(k)} = T_J^{\,k}\, e^{(0)}$ видно: метод сходится при любом $x^{(0)}$ тогда и только тогда, когда $T_J^{\,k} \to 0$ при $k \to \infty$. Это эквивалентно условию на спектральный радиус:

$$\rho(T_J) = \max_i |\lambda_i| < 1,$$

где $\lambda_i$ - собственные значения $T_J$, а $\rho$ - модуль наибольшего из них. Спектральный радиус - точный критерий, но его дорого считать. На практике пользуются достаточным условием - диагональным преобладанием по строкам:

$$|a_{ii}| > \sum_{j \neq i} |a_{ij}| \quad \text{для всех } i.$$

Если оно выполнено, то $\|T_J\|_\infty < 1$, спектральный радиус заведомо меньше единицы, и метод сходится при любом стартовом приближении. Условие лишь достаточное: метод может сходиться и без диагонального преобладания, но тогда сходимость нужно проверять через $\rho(T_J)$ напрямую. Подробнее общая теория простой итерации разобрана в статье про сходимость метода простой итерации.

Скорость сходимости

Ошибка убывает примерно как $\|e^{(k)}\| \approx \rho(T_J)^k \, \|e^{(0)}\|$. Чтобы уменьшить ошибку в $10$ раз, нужно около $1 / \log_{10}\!\big(1/\rho(T_J)\big)$ итераций. При $\rho = 0.5$ хватает $3$–$4$ шагов на порядок точности, а при $\rho = 0.99$ - уже около $230$. Именно поэтому Якоби пасует на задачах с плохо обусловленными матрицами: для дискретизации уравнения Пуассона на сетке $N \times N$ спектральный радиус $\rho(T_J) \approx 1 - O(h^2)$, где $h = 1/N$, и число итераций растёт как $O(N^2)$.

Для широкого класса согласованно упорядоченных матриц справедливо соотношение Янга $\rho(T_{GS}) = \rho(T_J)^2$, связывающее Якоби с Гауссом–Зейделем. Оно означает, что одна итерация Гаусса–Зейделя в среднем эквивалентна двум итерациям Якоби. Это главный практический минус Якоби: при одинаковой матрице он сходится примерно вдвое медленнее по числу шагов.

Сравнение с методом Гаусса–Зейделя

Метод Гаусса–Зейделя использует уже пересчитанные на текущем шаге компоненты $x_j^{(k+1)}$ при $j < i$, тогда как Якоби берёт только «старый» вектор. У этого есть две стороны.

• Скорость. Гаусс–Зейдель обычно сходится вдвое быстрее по числу итераций (соотношение Янга выше).
• Параллелизм. Якоби тривиально параллелится: все $n$ компонент независимы и считаются одновременно - идеально для GPU и векторных архитектур. Гаусс–Зейдель последователен «в лоб» и требует трюков вроде red–black упорядочения.
• Память. Якоби нужно хранить два вектора (старый $x^{(k)}$ и новый $x^{(k+1)}$), Гауссу–Зейделю - один, так как он переписывает массив на месте.

На многоядерных и графических процессорах выигрыш в параллелизме часто перевешивает проигрыш в числе итераций, поэтому Якоби (и его взвешенный вариант) до сих пор широко применяют как сглаживатель в многосеточных методах.

Критерий остановки

Точного решения мы не знаем, поэтому итерации прекращают, когда соседние приближения перестают заметно меняться. Корректный критерий - относительная норма поправки:

$$\frac{\|x^{(k+1)} - x^{(k)}\|}{\|x^{(k+1)}\|} < \varepsilon.$$

Часто используют норму невязки $\|b - A x^{(k)}\| < \varepsilon$ - она напрямую измеряет, насколько хорошо приближение удовлетворяет системе. Абсолютный критерий $\|x^{(k+1)} - x^{(k)}\| < \varepsilon$ без нормировки опасен: при больших по модулю компонентах решения он сработает слишком рано и даст ложную «сходимость».

Пошаговый пример

Решим систему $3 \times 3$:

$$\begin{cases} 10 x_1 - x_2 + 2 x_3 = 6, \\ -x_1 + 11 x_2 - x_3 = 25, \\ 2 x_1 - x_2 + 10 x_3 = -11. \end{cases}$$

Диагональное преобладание выполнено: $10 > |{-1}| + |2| = 3$, $11 > 2$, $10 > 3$ - метод сойдётся. Итерационные формулы:

$$x_1^{(k+1)} = \frac{6 + x_2^{(k)} - 2 x_3^{(k)}}{10}, \quad x_2^{(k+1)} = \frac{25 + x_1^{(k)} + x_3^{(k)}}{11}, \quad x_3^{(k+1)} = \frac{-11 - 2 x_1^{(k)} + x_2^{(k)}}{10}.$$

Стартуем с $x^{(0)} = (0, 0, 0)$. Первая итерация: $x_1^{(1)} = 0.6$, $x_2^{(1)} = 2.2727$, $x_3^{(1)} = -1.1$. Вторая: $x_1^{(2)} = 1.0473$, $x_2^{(2)} = 2.2273$, $x_3^{(2)} = -0.9927$. Значения уже близки к точному решению $x^{*} \approx (1.0433, 2.2692, -1.0817)$, и за $6$–$8$ итераций ошибка падает ниже $10^{-3}$.

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

• Применить метод к системе с нулём на диагонали без перестановки строк - формула $1/a_{ii}$ даст деление на ноль.
• Перепутать Якоби с Гауссом–Зейделем и подставлять уже обновлённые компоненты текущего шага - это меняет метод и его скорость сходимости.
• Объявить метод расходящимся после $5$–$10$ итераций: при $\rho(T_J)$ около $0.95$ требуются сотни шагов даже для умеренной точности.
• Считать, что диагональное преобладание необходимо. Это лишь достаточное условие - без него метод тоже может сходиться, нужно проверить $\rho(T_J) < 1$.
• Использовать абсолютный критерий остановки $\|x^{(k+1)} - x^{(k)}\| < \varepsilon$ без нормировки на больших по модулю решениях.

FAQ

Чем метод Якоби отличается от метода Гаусса–Зейделя? Якоби вычисляет все компоненты нового вектора $x^{(k+1)}$ только по «старому» вектору $x^{(k)}$, поэтому компоненты независимы и метод легко распараллеливается. Гаусс–Зейдель внутри одного шага сразу подставляет уже пересчитанные $x_j^{(k+1)}$ при $j < i$, за счёт чего сходится примерно вдвое быстрее, но теряет параллелизм.

Как проверить, сойдётся ли метод Якоби? Достаточное условие - диагональное преобладание по строкам: $|a_{ii}| > \sum_{j \neq i} |a_{ij}|$ для всех $i$. Если оно выполнено, метод сходится при любом $x^{(0)}$. Точный критерий - спектральный радиус матрицы перехода $\rho(T_J) < 1$, где $T_J = -D^{-1}(L + U)$.

Что делать, если диагональное преобладание не выполнено? Сначала попробовать переставить строки (и столбцы) так, чтобы наибольшие по модулю элементы оказались на диагонали. Если преобладания всё равно нет, посчитать спектральный радиус $\rho(T_J)$ численно: при $\rho < 1$ метод сойдётся, иначе нужен другой подход - Гаусс–Зейдель, SOR или сопряжённые градиенты.

Коротко

Метод Якоби - итерационный метод для $Ax = b$, где каждая компонента нового приближения считается из предыдущего вектора по формуле $x_i^{(k+1)} = (b_i - \sum_{j \neq i} a_{ij} x_j^{(k)}) / a_{ii}$. Сходимость гарантируется при $\rho(T_J) < 1$; удобное достаточное условие - диагональное преобладание. Метод проще и параллелится лучше Гаусса–Зейделя, но сходится примерно вдвое медленнее по числу итераций.