Метод Якоби

07.02.2021

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

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

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

A x → = b → {displaystyle A{vec {x}}={vec {b}}} , где A = ( a 11 … a 1 n ⋮ ⋱ ⋮ a n 1 … a n n ) , b → = ( b 1 ⋮ b n ) {displaystyle A=left({egin{array}{ccc}a_{11}&ldots &a_{1n}vdots &ddots &vdots a_{n1}&ldots &a_{nn}end{array}} ight),quad {vec {b}}=left({egin{array}{c}b_{1}vdots b_{n}end{array}} ight)}

Или { a 11 x 1 + … + a 1 n x n = b 1 …                                           a n 1 x 1 + … + a n n x n = b n {displaystyle left{{egin{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}} ight.}

Описание метода

Для того, чтобы построить итеративную процедуру метода Якоби, необходимо провести предварительное преобразование системы уравнений A x → = b → {displaystyle A{vec {x}}={vec {b}}} к итерационному виду x → = B x → + g → {displaystyle {vec {x}}=B{vec {x}}+{vec {g}}} . Оно может быть осуществлено по одному из следующих правил:

  • B = E − D − 1 A = D − 1 ( D − A ) , g → = D − 1 b → ; {displaystyle B=E-D^{-1}A=D^{-1}(D-A),quad {vec {g}}=D^{-1}{vec {b}};}
  • B = − D − 1 ( L + U ) = − D − 1 ( A − D ) , g → = D − 1 b → {displaystyle B=-D^{-1}(L+U)=-D^{-1}(A-D),quad {vec {g}}=D^{-1}{vec {b}}}
  • D i i − 1 = 1 / D i i , D i i ≠ 0 , i = 1 , 2 , . . . , n ; {displaystyle D_{ii}^{-1}=1/D_{ii},D_{ii} eq 0,,i=1,2,...,nquad ;}


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

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

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

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

x i ( k + 1 ) = 1 a i i ( b i − ∑ j ≠ i a i j x j ( k ) ) , i = 1 , 2 , … , n . {displaystyle x_{i}^{(k+1)}={frac {1}{a_{ii}}}left(b_{i}-sum _{j eq i}a_{ij}x_{j}^{(k)} ight),quad i=1,2,ldots ,n.}

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

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

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

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

Условие остановки

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

∥ x ( k + 1 ) − x ( k ) ∥≤ 1 − q q ε . {displaystyle parallel x^{(k+1)}-x^{(k)}parallel leq {frac {1-q}{q}}varepsilon .}

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

∥ x ( k + 1 ) − x ( k ) ∥≤ ε . {displaystyle parallel x^{(k+1)}-x^{(k)}parallel leq varepsilon .}

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

∥ A x ( k ) − b ∥≤ ε {displaystyle parallel Ax^{(k)}-bparallel 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; }