Метод Якоби

Возьмём систему линейных уравнений:

${\displaystyle A{\vec {x}}={\vec {b}}}$, где ${\displaystyle A=\left({\begin{array}{ccc}a_{11}&\ldots &a_{1n}\\\vdots &\ddots &\vdots \\a_{n1}&\ldots &a_{nn}\end{array}}\right),\quad {\vec {b}}=\left({\begin{array}{c}b_{1}\\\vdots \\b_{n}\end{array}}\right)}$

Или ${\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{\vec {x}}={\vec {b}}}$ к итерационному виду ${\displaystyle {\vec {x}}=B{\vec {x}}+{\vec {g}}}$. Оно может быть осуществлено по одному из следующих правил:

• ${\displaystyle B=E-D^{-1}A=D^{-1}(D-A),\quad {\vec {g}}=D^{-1}{\vec {b}};}$

• ${\displaystyle B=-D^{-1}(L+U)=-D^{-1}(A-D),\quad {\vec {g}}=D^{-1}{\vec {b}}}$

• ${\displaystyle D_{ii}^{-1}=1/D_{ii},D_{ii}\neq 0,\,i=1,2,...,n\quad$ ;}

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

Тогда процедура нахождения решения имеет вид:

${\displaystyle {\vec {x}}^{(k+1)}=B{\vec {x}}^{(k)}+{\vec {g}},}$

Или в виде поэлементной формулы:

${\displaystyle 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.}$

где ${\displaystyle k}$ счётчик итерации.

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

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

Теорема. Пусть ${\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 }$

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

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

#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;
}

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