Метод Якоби - Jacobi method

В числовой линейной алгебре , то метод Якоби представляет собой итерационный алгоритм для определения решений строго по диагонали доминирующей системы линейных уравнений . Решается для каждого диагонального элемента и подставляется приблизительное значение. Затем процесс повторяется до тех пор, пока он не сойдется. Этот алгоритм является урезанной версией метода диагонализации матриц с преобразованием Якоби . Метод назван в честь Карла Густава Якоба Якоби .

Описание

Позволять

${\ Displaystyle А \ mathbf {x} = \ mathbf {b}}$

- квадратная система из n линейных уравнений, где:

${\ displaystyle A = {\ begin {bmatrix} a_ {11} & a_ {12} & \ cdots & a_ {1n} \\ a_ {21} & a_ {22} & \ cdots & a_ {2n} \\\ vdots & \ vdots & \ ddots & \ vdots \\ a_ {n1} & a_ {n2} & \ cdots & a_ {nn} \ end {bmatrix}}, \ qquad \ mathbf {x} = {\ begin {bmatrix} x_ {1} \\ x_ {2} \\\ vdots \\ x_ {n} \ end {bmatrix}}, \ qquad \ mathbf {b} = {\ begin {bmatrix} b_ {1} \\ b_ {2} \\\ vdots \ \ b_ {n} \ end {bmatrix}}.}$

Тогда A можно разложить на диагональный компонент D , нижнюю треугольную часть L и верхнюю треугольную часть U :

${\ displaystyle A = D + L + U \ qquad {\ text {where}} \ qquad D = {\ begin {bmatrix} a_ {11} & 0 & \ cdots & 0 \\ 0 & a_ {22} & \ cdots & 0 \\\ vdots & \ vdots & \ ddots & \ vdots \\ 0 & 0 & \ cdots & a_ {nn} \ end {bmatrix}} {\ text {and}} L + U = {\ begin {bmatrix} 0 & a_ {12} & \ cdots & a_ {1n} \\ a_ {21} & 0 & \ cdots & a_ {2n} \\\ vdots & \ vdots & \ ddots & \ vdots \\ a_ {n1} & a_ {n2} & \ cdots & 0 \ end {bmatrix}}. }$

Затем решение итеративно получается через

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

где - k- е приближение или итерация, а - следующая или k + 1 итерация . Таким образом, формула на основе элементов: ${\ Displaystyle \ mathbf {х} ^ {(к)}}$${\ displaystyle \ mathbf {x}}$${\ Displaystyle \ mathbf {х} ^ {(к + 1)}}$${\ displaystyle \ mathbf {x}}$

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

Для вычисления требуется каждый элемент в x ^{( k ),} кроме самого себя. В отличие от метода Гаусса-Зейделя , мы не можем переписать с , так как это значение будет необходимо в остальной части вычисления. Минимальный объем памяти - два вектора размера n . ${\ Displaystyle х_ {я} ^ {(к + 1)}}$${\ Displaystyle х_ {я} ^ {(к)}}$${\ Displaystyle х_ {я} ^ {(к + 1)}}$

Алгоритм

Input: initial guess ${\displaystyle x^{(0)}}$ to the solution, (diagonal dominant) matrix ${\displaystyle A}$, right-hand side vector ${\displaystyle b}$, convergence criterion
Output: solution when convergence is reached
Comments: pseudocode based on the element-based formula above

${\displaystyle k=0}$
while convergence not reached do
    for i := 1 step until n do
        ${\displaystyle \sigma =0}$
        for j := 1 step until n do
            if j ≠ i then
                ${\displaystyle \sigma =\sigma +a_{ij}x_{j}^{(k)}}$
            end
        end
        ${\displaystyle x_{i}^{(k+1)}={{\frac {1}{a_{ii}}}\left({b_{i}-\sigma }\right)}}$
    end
    ${\displaystyle k=k+1}$
end

Конвергенция

Стандартное условие сходимости (для любого итерационного метода) - это когда спектральный радиус итерационной матрицы меньше 1:

${\ Displaystyle \ rho (D ^ {- 1} (L + U)) <1.}$

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

${\ displaystyle \ left | a_ {ii} \ right |> \ sum _ {j \ neq i} {\ left | a_ {ij} \ right |}.}$

Метод Якоби иногда сходится, даже если эти условия не выполняются.

Обратите внимание, что метод Якоби не сходится для каждой симметричной положительно определенной матрицы . Например

${\ displaystyle A = {\ begin {pmatrix} 29 & 2 & 1 \\ 2 & 6 & 1 \\ 1 & 1 & {\ frac {1} {5}} \ end {pmatrix}} \ quad \ Rightarrow \ quad D ^ {- 1} (L + U ) = {\ begin {pmatrix} 0 & {\ frac {2} {29}} & {\ frac {1} {29}} \\ {\ frac {1} {3}} & 0 & {\ frac {1} { 6}} \\ 5 & 5 & 0 \ end {pmatrix}} \ quad \ Rightarrow \ quad \ rho (D ^ {- 1} (L + U)) \ приблизительно 1.0661 \ ,.}$

Примеры

Пример 1

Линейная система вида с начальной оценкой дается формулой ${\ displaystyle Ax = b}$${\ displaystyle x ^ {(0)}}$

${\ displaystyle A = {\ begin {bmatrix} 2 & 1 \\ 5 & 7 \\\ end {bmatrix}}, \ b = {\ begin {bmatrix} 11 \\ 13 \\\ end {bmatrix}} \ quad {\ text {и}} \ quad x ^ {(0)} = {\ begin {bmatrix} 1 \\ 1 \\\ end {bmatrix}}.}$

Для оценки мы используем уравнение , описанное выше . Сначала перепишем уравнение в более удобном виде , где и . Из известных ценностей ${\ displaystyle x ^ {(k + 1)} = D ^ {- 1} (b- (L + U) x ^ {(k)})}$${\ displaystyle x}$${\ Displaystyle D ^ {- 1} (b- (L + U) x ^ {(k)}) = Tx ^ {(k)} + C}$${\ Displaystyle Т = -D ^ {- 1} (L + U)}$${\ displaystyle C = D ^ {- 1} b}$

${\ displaystyle D ^ {- 1} = {\ begin {bmatrix} 1/2 & 0 \\ 0 & 1/7 \\\ end {bmatrix}}, \ L = {\ begin {bmatrix} 0 & 0 \\ 5 & 0 \\\ конец {bmatrix}} \ quad {\ text {and}} \ quad U = {\ begin {bmatrix} 0 & 1 \\ 0 & 0 \\\ end {bmatrix}}.}$

мы определяем как ${\ Displaystyle Т = -D ^ {- 1} (L + U)}$

${\ displaystyle T = {\ begin {bmatrix} 1/2 & 0 \\ 0 & 1/7 \\\ end {bmatrix}} \ left \ {{\ begin {bmatrix} 0 & 0 \\ - 5 & 0 \\\ end {bmatrix}} + {\ begin {bmatrix} 0 & -1 \\ 0 & 0 \\\ end {bmatrix}} \ right \} = {\ begin {bmatrix} 0 & -1 / 2 \\ - 5/7 & 0 \\\ end {bmatrix} }.}$

Далее находится как ${\ displaystyle C}$

${\ displaystyle C = {\ begin {bmatrix} 1/2 & 0 \\ 0 & 1/7 \\\ end {bmatrix}} {\ begin {bmatrix} 11 \\ 13 \\\ end {bmatrix}} = {\ begin { bmatrix} 11/2 \\ 13/7 \\\ end {bmatrix}}.}$

Имея и рассчитав, мы оцениваем как : ${\ displaystyle T}$${\ displaystyle C}$${\ displaystyle x}$${\ Displaystyle х ^ {(1)} = Tx ^ {(0)} + C}$

${\ displaystyle x ^ {(1)} = {\ begin {bmatrix} 0 & -1/2 \\ - 5/7 & 0 \\\ end {bmatrix}} {\ begin {bmatrix} 1 \\ 1 \\\ конец {bmatrix}} + {\ begin {bmatrix} 11/2 \\ 13/7 \\\ end {bmatrix}} = {\ begin {bmatrix} 5.0 \\ 8/7 \\\ end {bmatrix}} \ приблизительно {\ begin {bmatrix} 5 \\ 1.143 \\\ end {bmatrix}}.}$

Следующая итерация дает

${\ displaystyle x ^ {(2)} = {\ begin {bmatrix} 0 & -1/2 \\ - 5/7 & 0 \\\ end {bmatrix}} {\ begin {bmatrix} 5.0 \\ 8/7 \\ \ end {bmatrix}} + {\ begin {bmatrix} 11/2 \\ 13/7 \\\ end {bmatrix}} = {\ begin {bmatrix} 69/14 \\ - 12/7 \\\ end { bmatrix}} \ приблизительно {\ begin {bmatrix} 4,929 \\ - 1,714 \\\ end {bmatrix}}.}$

Этот процесс повторяется до схождения (т. Е. До тех пор, пока не станет маленьким). Решение после 25 итераций: ${\ displaystyle \ | Ax ^ {(n)} - b \ |}$

${\ displaystyle x = {\ begin {bmatrix} 7.111 \\ - 3.222 \ end {bmatrix}}.}$

Пример 2

Предположим, нам дана следующая линейная система:

${\ displaystyle {\ begin {align} 10x_ {1} -x_ {2} + 2x_ {3} & = 6, \\ - x_ {1} + 11x_ {2} -x_ {3} + 3x_ {4} & = 25, \\ 2x_ {1} -x_ {2} + 10x_ {3} -x_ {4} & = - 11, \\ 3x_ {2} -x_ {3} + 8x_ {4} & = 15. \ конец {выровнен}}}$

Если мы выберем (0, 0, 0, 0) в качестве начального приближения, то первое приближенное решение будет иметь вид

${\ displaystyle {\ begin {align} x_ {1} & = (6 + 0- (2 * 0)) / 10 = 0,6, \\ x_ {2} & = (25 + 0 + 0- (3 * 0 )) / 11 = 25/11 = 2,2727, \\ x_ {3} & = (- 11- (2 * 0) + 0 + 0) /10=-1,1, \\ x_ {4} & = (15- (3 * 0) +0) /8=1.875.\end {выровнено}}}$

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

${\ displaystyle x_ {1}}$	${\ displaystyle x_ {2}}$	${\ displaystyle x_ {3}}$	${\ displaystyle x_ {4}}$
0,6	2,27272	-1,1	1,875
1,04727	1,7159	-0,80522	0,88522
0,93263	2,05330	-1,0493	1,13088
1.01519	1,95369	-0,9681	0,97384
0,98899	2,0114	-1,0102	1,02135

Точное решение системы - (1, 2, −1, 1) .

Пример Python

import numpy as np

ITERATION_LIMIT = 1000

# initialize the matrix
A = np.array([[10., -1., 2., 0.],
              [-1., 11., -1., 3.],
              [2., -1., 10., -1.],
              [0.0, 3., -1., 8.]])
# initialize the RHS vector
b = np.array([6., 25., -11., 15.])

# prints the system
print("System:")
for i in range(A.shape[0]):
    row = ["{}*x{}".format(A[i, j], j + 1) for j in range(A.shape[1])]
    print(" + ".join(row), "=", b[i])
print()

x = np.zeros_like(b)
for it_count in range(ITERATION_LIMIT):
    if it_count != 0:
        print("Iteration {0}: {1}".format(it_count, x))
    x_new = np.zeros_like(x)

    for i in range(A.shape[0]):
        s1 = np.dot(A[i, :i], x[:i])
        s2 = np.dot(A[i, i + 1:], x[i + 1:])
        x_new[i] = (b[i] - s1 - s2) / A[i, i]
        if x_new[i] == x_new[i-1]:
          break

    if np.allclose(x, x_new, atol=1e-10, rtol=0.):
        break

    x = x_new

print("Solution:")
print(x)
error = np.dot(A, x) - b
print("Error:")
print(error)

Весовой метод Якоби

Взвешенная итерация Якоби использует параметр для вычисления итерации как ${\ displaystyle \ omega}$

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

с обычным выбором. Из отношения это также может быть выражено как ${\ displaystyle \ omega = 2/3}$${\ Displaystyle L + U = AD}$

${\ displaystyle \ mathbf {x} ^ {(k + 1)} = \ omega D ^ {- 1} \ mathbf {b} + \ left (I- \ omega D ^ {- 1} A \ right) \ mathbf {x} ^ {(k)}}$.

Сходимость в симметричном положительно определенном случае

В случае, если матрица системы имеет симметричный положительно определенный тип, можно показать сходимость. ${\ displaystyle A}$

Позвольте быть итерационной матрицей. Тогда сходимость гарантирована для ${\ displaystyle C = C _ {\ omega} = I- \ omega D ^ {- 1} A}$

${\ displaystyle \ rho (C _ {\ omega}) <1 \ quad \ Longleftrightarrow \ quad 0 <\ omega <{\ frac {2} {\ lambda _ {\ text {max}} (D ^ {- 1} A )}} \ ,,}$

где - максимальное собственное значение. ${\ displaystyle \ lambda _ {\ text {max}}}$

Спектральный радиус может быть минимизирован для конкретного выбора следующим образом: ${\ displaystyle \ omega = \ omega _ {\ text {opt}}}$

${\ displaystyle \ min _ {\ omega} \ rho (C _ {\ omega}) = \ rho (C _ {\ omega _ {\ text {opt}}}) = 1 - {\ frac {2} {\ kappa ( D ^ {- 1} A) +1}} \ quad {\ text {for}} \ quad \ omega _ {\ text {opt}}: = {\ frac {2} {\ lambda _ {\ text {min }} (D ^ {- 1} A) + \ lambda _ {\ text {max}} (D ^ {- 1} A)}} \ ,,}$

где - номер обусловленности матрицы . ${\ displaystyle \ kappa}$

Смотрите также

• Метод Гаусса – Зейделя.
• Последовательное чрезмерное расслабление
• Итерационный метод § Линейные системы
• Распространение гауссовской веры
• Расщепление матрицы

использованная литература

внешние ссылки

• Эта статья включает текст из статьи Jacobi_method на CFD-Wiki, которая находится под лицензией GFDL .

• Блэк, Ноэль; Мур, Ширли и Вайстейн, Эрик В. «Метод Якоби» . MathWorld .
• Метод Якоби с сайта www.math-linux.com