Метод Якоби

Метод Якоби — разновидность метода простой итерации для численного решения системы линейных алгебраических уравнений. Назван в честь Карла Густава Якоби.

Пусть требуется численно решить систему линейных уравнений:

${\displaystyle \left\{{\begin{array}{rcl}a_{11}x_{1}+\ldots +a_{1n}x_{n}&=&b_{1}\\\ \ldots ~~~~~~~~~~~~~~~~~~~~~~~\\a_{n1}x_{1}+\ldots +a_{nn}x_{n}&=&b_{n}\end{array}}\right.}$

Предполагается, что ${\displaystyle a_{ii}\neq 0,\ i\in \{1,\ldots ,n\}}$ (иначе метод Якоби неприменим). Выразим ${\displaystyle x_{1}}$ через первое уравнение, ${\displaystyle x_{2}}$ — через второе и т. д.^[1]:

${\displaystyle \left\{{\begin{array}{rcl}x_{1}=&{\dfrac {1}{a_{11}}}\left(b_{1}-a_{12}x_{2}-\ldots -a_{1n}x_{n}\right)\\\ldots \\x_{n}=&{\dfrac {1}{a_{nn}}}\left(b_{n}-a_{n1}x_{1}-\ldots -a_{n(n-1)}x_{n-1}\right)\end{array}}\right.}$

В методе Якоби последовательность приближений ${\displaystyle x^{(k)}}$ строится следующим образом. Выбирается первое приближение ${\displaystyle x^{(0)}}$, формула для остальных приближений имеет вид

${\displaystyle {\begin{array}{rcl}x_{1}^{(k+1)}=&{\dfrac {1}{a_{11}}}\left(b_{1}-a_{12}x_{2}^{(k)}-\ldots -a_{1n}x_{n}^{(k)}\right)\\\ldots \\x_{n}^{(k+1)}=&{\dfrac {1}{a_{nn}}}\left(b_{n}-a_{n1}x_{1}^{(k)}-\ldots -a_{n(n-1)}x_{n-1}^{(k)}\right)\end{array}}}$.

В матричной форме имеет следующий вид. Пусть СЛАУ в матричной форме записано как

${\displaystyle Ax=b}$

Представим матрицу ${\displaystyle A}$ в виде ${\displaystyle A=L+D+U}$, где ${\displaystyle D}$ означает диагональную матрицу, у которой на главной диагонали стоят соответствующие элементы матрицы ${\displaystyle A}$; тогда как матрицы ${\displaystyle U}$ и ${\displaystyle L}$ содержат верхнюю и нижнюю треугольные части ${\displaystyle A}$, на главной диагонали которых нули. Тогда итерационную формулу можно записать как

${\displaystyle x^{(k+1)}=D^{-1}(b-(L+U)x^{(k)})}$

Приведем достаточное условие сходимости метода.

Теорема. Пусть ${\displaystyle \|B\|<1}$. Тогда при любом выборе начального приближения ${\displaystyle {\vec {x}}^{(0)}}$: метод сходится; скорость сходимости метода равна скорости сходимости геометрической прогрессии со знаменателем ${\displaystyle q=\|B\|}$; верна оценка погрешности: ${\displaystyle \|{\vec {x}}^{(k)}-{\vec {x}}\|<q^{k}\,\|{\vec {x}}^{(0)}-{\vec {x}}\|}$.

Условие окончания итерационного процесса при достижении точности ${\displaystyle \varepsilon }$ имеет вид:

${\displaystyle \parallel x^{(k+1)}-x^{(k)}\parallel \leq {\frac {1-q}{q}}\varepsilon .}$

Для достаточно хорошо обусловленной матрицы ${\displaystyle B}$ (при ${\displaystyle \parallel B\parallel =q<1/2}$) оно выполняется при

${\displaystyle \parallel x^{(k+1)}-x^{(k)}\parallel \leq \varepsilon .}$

Данный критерий не требует вычисления нормы матрицы и часто используется на практике. При этом точное условие окончания итерационного процесса имеет вид

${\displaystyle \parallel Ax^{(k)}-b\parallel \leq \varepsilon }$

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

В отличие от метода Гаусса-Зейделя мы не можем заменять ${\displaystyle x_{i}^{(k)}}$ на ${\displaystyle x_{i}^{(k+1)}}$ в процессе итерационной процедуры, так как эти значения понадобятся для остальных вычислений. Это наиболее значимое различие между методом Якоби и методом Гаусса-Зейделя решения СЛАУ. Таким образом на каждой итерации придётся хранить оба вектора приближений: старый и новый.

Ниже приведён алгоритм реализации на C++

#include <math.h> const double eps = 0.001; ///< желаемая точность   ..........................  /// N - размерность матрицы; A[N][N] - матрица коэффициентов, F[N] - столбец свободных членов, /// X[N] - начальное приближение, ответ записывается также в X[N]; void Jacobi (int N, double** A, double* F, double* X) { 	double* TempX = new double[N]; 	double norm; // норма, определяемая как наибольшая разность компонент столбца иксов соседних итераций.  	do { 		for (int i = 0; i < N; i++) { 			TempX[i] = F[i]; 			for (int g = 0; g < N; g++) { 				if (i != g) 					TempX[i] -= A[i][g] * X[g]; 			} 			TempX[i] /= A[i][i]; 		}         norm = fabs(X[0] - TempX[0]); 		for (int h = 0; h < N; h++) { 			if (fabs(X[h] - TempX[h]) > norm) 				norm = fabs(X[h] - TempX[h]); 			X[h] = TempX[h]; 		} 	} while (norm > eps); 	delete[] TempX; } 

Ниже приведен алгоритм реализации на Python

from collections.abc import Sequence, MutableSequence   def Jacobi(         A: Sequence[Sequence[float]],         b: Sequence[float],         eps: float = 0.001,         x_init: MutableSequence[float] | None = None) -> list[float]:     """     метод Якоби для A*x = b (*)      :param A: Матрица (*) слева      :param b: известный вектор (*) справа      :param x_init: начальное приближение      :return: приблизительное значения х (*)     """      def sum(a: Sequence[float], x: Sequence[float], j: int) -> float:         S: float = 0         for i, (m, y) in enumerate(zip(a, x)):             if i != j:                 S += m*y         return S      def norm(x: Sequence[float], y: Sequence[float]) -> float:         return max((abs(x0-y0) for x0, y0 in zip(x, y)))      if x_init is None:         y = [a/A[i][i] for i, a in enumerate(b)]     else:         y = x.copy()      x: list[float] = [-(sum(a, y, i) - b[i])/A[i][i]                       for i, a in enumerate(A)]      while norm(y, x) > eps:         for i, elem in enumerate(x):             y[i] = elem         for i, a in enumerate(A):             x[i] = -(sum(a, y, i) - b[i])/A[i][i]     return x 

• Березин, И. С., Жидков Н. П. Методы вычислений (рус.). — М.: Физматлит, 1959. — Т. II.

• Метод
• Метод Гаусса — Зейделя

Для улучшения этой статьи по математике желательно: Найти и оформить в виде сносок ссылки на независимые авторитетные источники, подтверждающие написанное. После исправления проблемы исключите её из списка. Удалите шаблон, если устранены все недостатки.