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

В числовой линейной алгебре метод Гаусса –Зейделя , также известный как метод Либмана или метод последовательных смещений , представляет собой итерационный метод, используемый для решения системы линейных уравнений . Он назван в честь немецких математиков Карла Фридриха Гаусса и Филиппа Людвига фон Зейделя . Хотя его можно применить к любой матрице с ненулевыми элементами на диагоналях, сходимость гарантируется только в том случае, если матрица либо строго доминирует по диагонали , или симметричный и положительно определенный . Об этом упоминается только в частном письме Гаусса своему ученику Герлингу в 1823 году. Публикация была выпущена Зейделем только в 1874 году.

Описание

Позволять ${\textstyle \mathbf {A} \mathbf {x} =\mathbf {b} }$ — квадратная система n линейных уравнений, где:

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

Когда $\mathbf {A}$ и $\mathbf {b}$ известны, и $\mathbf {x}$ неизвестно, мы можем использовать метод Гаусса – Зейделя для аппроксимации $\mathbf {x}$. Вектор $\mathbf {x} ^{(0)}$ обозначает наше первоначальное предположение для $\mathbf {x}$ (часто $\mathbf {x} _{i}^{(0)}=0$ для $i=1,2,...,n$). Обозначим через $\mathbf {x} ^{(k)}$ тот $k$-е приближение или итерация $\mathbf {x}$и по $\mathbf {x} ^{(k+1)}$ приближение $\mathbf {x}$ на следующем (или $k+1$) итерация.

Формула на основе матрицы

Решение получается итерационно через $$\mathbf {L} \mathbf {x} ^{(k+1)}=\mathbf {b} -\mathbf {U} \mathbf {x} ^{(k)},$$ где матрица $\mathbf {A}$ разлагается на нижнюю треугольную составляющую $\mathbf {L}$и строго верхнетреугольная компонента $\mathbf {U}$ такой, что $\mathbf {A} =\mathbf {L} +\mathbf {U}$. Точнее, разложение $A$ в $L_{*}$ и $U$ дается:

$$\mathbf {A} =\underbrace {\begin{bmatrix}a_{11}&0&\cdots &0\\a_{21}&a_{22}&\cdots &0\\\vdots &\vdots &\ddots &\vdots \\a_{n1}&a_{n2}&\cdots &a_{nn}\end{bmatrix}} _{\textstyle \mathbf {L} }+\underbrace {\begin{bmatrix}0&a_{12}&\cdots &a_{1n}\\0&0&\cdots &a_{2n}\\\vdots &\vdots &\ddots &\vdots \\0&0&\cdots &0\end{bmatrix}} _{\textstyle \mathbf {U} }.$$

Почему работает матричная формула

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

${\begin{alignedat}{1}\mathbf {A} \mathbf {x} &=\mathbf {b} \\(\mathbf {L} +\mathbf {U} )\mathbf {x} &=\mathbf {b} \\\mathbf {L} \mathbf {x} +\mathbf {U} \mathbf {x} &=\mathbf {b} \\\mathbf {L} \mathbf {x} &=\mathbf {b} -\mathbf {U} \mathbf {x} \end{alignedat}}$

Метод Гаусса – Зейделя теперь решает левую часть этого выражения для $\mathbf {x}$, используя предыдущее значение для $\mathbf {x}$ с правой стороны. Аналитически это можно записать как $$\mathbf {x} ^{(k+1)}=\mathbf {L} ^{-1}\left(\mathbf {b} -\mathbf {U} \mathbf {x} ^{(k)}\right).$$

Формула на основе элементов

Однако, воспользовавшись треугольной формой $\mathbf {L}$, элементы $\mathbf {x} ^{(k+1)}$ можно вычислить последовательно для каждой строки $i$ используя прямую замену : $$x_{i}^{(k+1)}={\frac {1}{a_{ii}}}\left(b_{i}-\sum _{j=1}^{i-1}a_{ij}x_{j}^{(k+1)}-\sum _{j=i+1}^{n}a_{ij}x_{j}^{(k)}\right),\quad i=1,2,\dots ,n.$$

Обратите внимание, что в формуле используются два суммирования на итерацию, которые можно выразить как одно суммирование. $\sum _{j\neq i}a_{ij}x_{j}$ который использует самую последнюю вычисленную итерацию $x_{j}$. Процедура обычно продолжается до тех пор, пока изменения, внесенные в ходе итерации, не окажутся ниже некоторого допуска, например, достаточно малого остатка .

Обсуждение

Поэлементная формула метода Гаусса – Зейделя связана с формулой (итеративного) метода Якоби с важным отличием:

В методе Гаусса-Зейделя вычисление $\mathbf {x} ^{(k+1)}$ использует элементы $\mathbf {x} ^{(k+1)}$ которые уже вычислены, и только элементы $\mathbf {x} ^{(k)}$ которые не были рассчитаны в $(k+1)$-я итерация. Это означает, что, в отличие от метода Якоби, требуется только один вектор хранения, поскольку элементы могут перезаписываться по мере их вычисления, что может быть полезно для очень больших задач.

Однако, в отличие от метода Якоби, вычисления для каждого элемента, как правило, гораздо сложнее реализовать параллельно , поскольку они могут иметь очень длинный критический путь и, следовательно, наиболее осуществимы для разреженных матриц . Более того, значения на каждой итерации зависят от порядка исходных уравнений.

Гаусса-Зейделя — это то же самое, что и последовательное чрезмерное расслабление с $\omega =1$.

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

Свойства сходимости метода Гаусса–Зейделя зависят от матрицы $\mathbf {A}$. А именно, известно, что процедура сходится, если:

• $\mathbf {A}$ является симметричным положительно определенным , или
• $\mathbf {A}$ строго или неприводимо доминирует по диагонали .

Метод Гаусса–Зейделя иногда сходится, даже если эти условия не выполняются.

Голуб и Ван Лоан приводят теорему для алгоритма, который расщепляет $\mathbf {A}$ на две части. Предполагать $\mathbf {A} =\mathbf {M} -\mathbf {N}$ является неособым. Позволять $r=\rho (\mathbf {M} ^{-1}\mathbf {N} )$ быть радиусом спектральным $\mathbf {M} ^{-1}\mathbf {N}$. Затем повторяется $\mathbf {x} ^{(k)}$ определяется $\mathbf {M} \mathbf {x} ^{(k+1)}=\mathbf {N} \mathbf {x} ^{(k)}+\mathbf {b}$ сходиться к $\mathbf {x} =\mathbf {A} ^{-1}\mathbf {b}$ для любого стартового вектора $\mathbf {x} ^{(0)}$ если $\mathbf {M}$ является неособым и $r<1$.

Алгоритм

Поскольку элементы могут быть перезаписаны во время их вычисления в этом алгоритме, требуется только один вектор хранения, а индексирование вектора опускается. Алгоритм следующий:

algorithm Gauss–Seidel method is
    inputs: A, b
    output: φ

    Choose an initial guess φ to the solution
    repeat until convergence
        for i from 1 until n do
            σ ← 0
            for j from 1 until n do
                if j ≠ i then
                    σ ← σ + aijφj
                end if
            end (j-loop)
            φi ← (bi − σ) / aii
        end (i-loop)
        check if convergence is reached
    end (repeat)

Примеры

Пример для матричной версии

Линейная система, показанная как $\mathbf {A} \mathbf {x} =\mathbf {b}$ дается: $$\mathbf {A} ={\begin{bmatrix}16&3\\7&-11\\\end{bmatrix}}\quad {\text{and}}\quad \mathbf {b} ={\begin{bmatrix}11\\13\end{bmatrix}}.$$

Мы хотим использовать уравнение $$\mathbf {x} ^{(k+1)}=\mathbf {L} ^{-1}(\mathbf {b} -\mathbf {U} \mathbf {x} ^{(k)})$$ в форме $$\mathbf {x} ^{(k+1)}=\mathbf {T} \mathbf {x} ^{(k)}+\mathbf {c}$$ где:

$\mathbf {T} =-\mathbf {L} ^{-1}\mathbf {U} \quad {\text{and}}\quad \mathbf {c} =\mathbf {L} ^{-1}\mathbf {b} .$

Мы должны разложить $\mathbf {A}$ в сумму нижней треугольной составляющей $\mathbf {L}$ и строгая верхняя треугольная составляющая $U$: $$\mathbf {L} ={\begin{bmatrix}16&0\\7&-11\\\end{bmatrix}}\quad {\text{and}}\quad \mathbf {U} ={\begin{bmatrix}0&3\\0&0\end{bmatrix}}.$$

Обратная сторона $\mathbf {L}$ является: $$\mathbf {L} ^{-1}={\begin{bmatrix}16&0\\7&-11\end{bmatrix}}^{-1}={\begin{bmatrix}0.0625&0.0000\\0.0398&-0.0909\\\end{bmatrix}}.$$

Теперь мы можем найти: $${\begin{aligned}\mathbf {T} &=-{\begin{bmatrix}0.0625&0.0000\\0.0398&-0.0909\end{bmatrix}}{\begin{bmatrix}0&3\\0&0\end{bmatrix}}={\begin{bmatrix}0.000&-0.1875\\0.000&-0.1194\end{bmatrix}},\\[1ex]\mathbf {c} &={\begin{bmatrix}0.0625&0.0000\\0.0398&-0.0909\end{bmatrix}}{\begin{bmatrix}11\\13\end{bmatrix}}={\begin{bmatrix}0.6875\\-0.7439\end{bmatrix}}.\end{aligned}}$$

Теперь у нас есть $\mathbf {T}$ и $\mathbf {c}$ и мы можем использовать их для получения векторов $\mathbf {x}$ итеративно.

Прежде всего, нам предстоит выбрать $\mathbf {x} ^{(0)}$: мы можем только догадываться. Чем точнее предположение, тем быстрее будет работать алгоритм.

Выбираем отправную точку: $$\mathbf {x} ^{(0)}={\begin{bmatrix}1.0\\1.0\end{bmatrix}}.$$

Затем мы можем рассчитать: $${\begin{aligned}\mathbf {x} ^{(1)}&={\begin{bmatrix}0.000&-0.1875\\0.000&-0.1193\end{bmatrix}}{\begin{bmatrix}1.0\\1.0\end{bmatrix}}+{\begin{bmatrix}0.6875\\-0.7443\end{bmatrix}}={\begin{bmatrix}0.5000\\-0.8636\end{bmatrix}}.\\[1ex]\mathbf {x} ^{(2)}&={\begin{bmatrix}0.000&-0.1875\\0.000&-0.1193\end{bmatrix}}{\begin{bmatrix}0.5000\\-0.8636\end{bmatrix}}+{\begin{bmatrix}0.6875\\-0.7443\end{bmatrix}}={\begin{bmatrix}0.8494\\-0.6413\end{bmatrix}}.\\[1ex]\mathbf {x} ^{(3)}&={\begin{bmatrix}0.000&-0.1875\\0.000&-0.1193\end{bmatrix}}{\begin{bmatrix}0.8494\\-0.6413\\\end{bmatrix}}+{\begin{bmatrix}0.6875\\-0.7443\end{bmatrix}}={\begin{bmatrix}0.8077\\-0.6678\end{bmatrix}}.\\[1ex]\mathbf {x} ^{(4)}&={\begin{bmatrix}0.000&-0.1875\\0.000&-0.1193\end{bmatrix}}{\begin{bmatrix}0.8077\\-0.6678\end{bmatrix}}+{\begin{bmatrix}0.6875\\-0.7443\end{bmatrix}}={\begin{bmatrix}0.8127\\-0.6646\end{bmatrix}}.\\[1ex]\mathbf {x} ^{(5)}&={\begin{bmatrix}0.000&-0.1875\\0.000&-0.1193\end{bmatrix}}{\begin{bmatrix}0.8127\\-0.6646\end{bmatrix}}+{\begin{bmatrix}0.6875\\-0.7443\end{bmatrix}}={\begin{bmatrix}0.8121\\-0.6650\end{bmatrix}}.\\[1ex]\mathbf {x} ^{(6)}&={\begin{bmatrix}0.000&-0.1875\\0.000&-0.1193\end{bmatrix}}{\begin{bmatrix}0.8121\\-0.6650\end{bmatrix}}+{\begin{bmatrix}0.6875\\-0.7443\end{bmatrix}}={\begin{bmatrix}0.8122\\-0.6650\end{bmatrix}}.\\[1ex]\mathbf {x} ^{(7)}&={\begin{bmatrix}0.000&-0.1875\\0.000&-0.1193\end{bmatrix}}{\begin{bmatrix}0.8122\\-0.6650\end{bmatrix}}+{\begin{bmatrix}0.6875\\-0.7443\end{bmatrix}}={\begin{bmatrix}0.8122\\-0.6650\end{bmatrix}}.\end{aligned}}$$

Как и ожидалось, алгоритм сходится к решению: $$\mathbf {x} =\mathbf {A} ^{-1}\mathbf {b} \approx {\begin{bmatrix}0.8122\\-0.6650\end{bmatrix}}.$$

Фактически, матрица A является строго диагонально-доминантной (но не положительно определенной).

Еще один пример для матричной версии

Другая линейная система, показанная как $\mathbf {A} \mathbf {x} =\mathbf {b}$ дается:

$$\mathbf {A} ={\begin{bmatrix}2&3\\5&7\\\end{bmatrix}}\quad {\text{and}}\quad \mathbf {b} ={\begin{bmatrix}11\\13\\\end{bmatrix}}.$$

Мы хотим использовать уравнение $$\mathbf {x} ^{(k+1)}=\mathbf {L} ^{-1}(\mathbf {b} -\mathbf {U} \mathbf {x} ^{(k)})$$ в форме $$\mathbf {x} ^{(k+1)}=\mathbf {T} \mathbf {x} ^{(k)}+\mathbf {c}$$ где:

$\mathbf {T} =-\mathbf {L} ^{-1}\mathbf {U} \quad {\text{and}}\quad \mathbf {c} =\mathbf {L} ^{-1}\mathbf {b} .$

Мы должны разложить $\mathbf {A}$ в сумму нижней треугольной составляющей $\mathbf {L}$ и строгая верхняя треугольная составляющая $\mathbf {U}$: $$\mathbf {L} ={\begin{bmatrix}2&0\\5&7\\\end{bmatrix}}\quad {\text{and}}\quad \mathbf {U} ={\begin{bmatrix}0&3\\0&0\\\end{bmatrix}}.$$

Обратная сторона $\mathbf {L}$ является: $$\mathbf {L} ^{-1}={\begin{bmatrix}2&0\\5&7\\\end{bmatrix}}^{-1}={\begin{bmatrix}0.500&0.000\\-0.357&0.143\\\end{bmatrix}}.$$

Теперь мы можем найти: $${\begin{aligned}\mathbf {T} &=-{\begin{bmatrix}0.500&0.000\\-0.357&0.143\\\end{bmatrix}}{\begin{bmatrix}0&3\\0&0\\\end{bmatrix}}={\begin{bmatrix}0.000&-1.500\\0.000&1.071\\\end{bmatrix}},\\[1ex]\mathbf {c} &={\begin{bmatrix}0.500&0.000\\-0.357&0.143\\\end{bmatrix}}{\begin{bmatrix}11\\13\\\end{bmatrix}}={\begin{bmatrix}5.500\\-2.071\\\end{bmatrix}}.\end{aligned}}$$

Теперь у нас есть $\mathbf {T}$ и $\mathbf {c}$ и мы можем использовать их для получения векторов $\mathbf {x}$ итеративно.

Прежде всего, нам предстоит выбрать $\mathbf {x} ^{(0)}$: мы можем только догадываться. Чем точнее предположение, тем быстрее будет работать алгоритм.

Мы предполагаем: $$\mathbf {x} ^{(0)}={\begin{bmatrix}1.1\\2.3\end{bmatrix}}.$$

Затем мы можем рассчитать: $${\begin{aligned}\mathbf {x} ^{(1)}&={\begin{bmatrix}0&-1.500\\0&1.071\\\end{bmatrix}}{\begin{bmatrix}1.1\\2.3\\\end{bmatrix}}+{\begin{bmatrix}5.500\\-2.071\\\end{bmatrix}}={\begin{bmatrix}2.050\\0.393\\\end{bmatrix}}.\\[1ex]\mathbf {x} ^{(2)}&={\begin{bmatrix}0&-1.500\\0&1.071\\\end{bmatrix}}{\begin{bmatrix}2.050\\0.393\\\end{bmatrix}}+{\begin{bmatrix}5.500\\-2.071\\\end{bmatrix}}={\begin{bmatrix}4.911\\-1.651\end{bmatrix}}.\\[1ex]\mathbf {x} ^{(3)}&=\cdots .\end{aligned}}$$

Если мы проверим сходимость, мы обнаружим, что алгоритм расходится. Фактически, матрица $\mathbf {A}$ не является ни диагонально-доминантным, ни положительно-определенным. Тогда сходимость к точному решению $$\mathbf {x} =\mathbf {A} ^{-1}\mathbf {b} ={\begin{bmatrix}-38\\29\end{bmatrix}}$$ не гарантировано и в данном случае не произойдет.

Пример версии уравнения

Предположим, что дано $n$ уравнения и отправная точка $\mathbf {x} _{0}$. На любом этапе итерации Гаусса-Зейделя решите первое уравнение для $x_{1}$ с точки зрения $x_{2},\dots ,x_{n}$; затем решите второе уравнение для $x_{2}$ с точки зрения $x_{1}$ только что нашел и остальное $x_{3},\dots ,x_{n}$; и продолжать $x_{n}$. Затем повторяйте итерации до тех пор, пока (надеюсь) не сойдемся.

Чтобы было понятно, рассмотрим пример: $${\begin{array}{rrrrl}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.\end{array}}$$

Решение для $x_{1},x_{2},x_{3}$ и $x_{4}$ дает: $${\begin{aligned}x_{1}&=x_{2}/10-x_{3}/5+3/5,\\x_{2}&=x_{1}/11+x_{3}/11-3x_{4}/11+25/11,\\x_{3}&=-x_{1}/5+x_{2}/10+x_{4}/10-11/10,\\x_{4}&=-3x_{2}/8+x_{3}/8+15/8.\end{aligned}}$$

Предположим, мы выбрали (0, 0, 0, 0) в качестве начального приближения, тогда первое приближенное решение имеет вид $${\begin{aligned}x_{1}&=3/5=0.6,\\x_{2}&=(3/5)/11+25/11=3/55+25/11=2.3272,\\x_{3}&=-(3/5)/5+(2.3272)/10-11/10=-3/25+0.23272-1.1=-0.9873,\\x_{4}&=-3(2.3272)/8+(-0.9873)/8+15/8=0.8789.\end{aligned}}$$

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

$x_{1}$	$x_{2}$	$x_{3}$	$x_{4}$
0.6	2.32727	−0.987273	0.878864
1.03018	2.03694	−1.01446	0.984341
1.00659	2.00356	−1.00253	0.998351
1.00086	2.0003	−1.00031	0.99985

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

Пример использования Python и NumPy

Следующая численная процедура просто повторяется для получения вектора решения.

import numpy as np

ITERATION_LIMIT = 1000

# initialize the matrix
A = np.array(
    [
        [10.0, -1.0, 2.0, 0.0],
        [-1.0, 11.0, -1.0, 3.0],
        [2.0, -1.0, 10.0, -1.0],
        [0.0, 3.0, -1.0, 8.0],
    ]
)
# initialize the RHS vector
b = np.array([6.0, 25.0, -11.0, 15.0])

print("System of equations:")
for i in range(A.shape[0]):
    row = [f"{A[i,j]:3g}*x{j+1}" for j in range(A.shape[1])]
    print("[{0}] = [{1:3g}]".format(" + ".join(row), b[i]))

x = np.zeros_like(b, np.float_)
for it_count in range(1, ITERATION_LIMIT):
    x_new = np.zeros_like(x, dtype=np.float_)
    print(f"Iteration {it_count}: {x}")
    for i in range(A.shape[0]):
        s1 = np.dot(A[i, :i], x_new[:i])
        s2 = np.dot(A[i, i + 1 :], x[i + 1 :])
        x_new[i] = (b[i] - s1 - s2) / A[i, i]
    if np.allclose(x, x_new, rtol=1e-8):
        break
    x = x_new

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

Производит вывод:

System of equations:
[ 10*x1 +  -1*x2 +   2*x3 +   0*x4] = [  6]
[ -1*x1 +  11*x2 +  -1*x3 +   3*x4] = [ 25]
[  2*x1 +  -1*x2 +  10*x3 +  -1*x4] = [-11]
[  0*x1 +   3*x2 +  -1*x3 +   8*x4] = [ 15]
Iteration 1: [ 0.  0.  0.  0.]
Iteration 2: [ 0.6         2.32727273 -0.98727273  0.87886364]
Iteration 3: [ 1.03018182  2.03693802 -1.0144562   0.98434122]
Iteration 4: [ 1.00658504  2.00355502 -1.00252738  0.99835095]
Iteration 5: [ 1.00086098  2.00029825 -1.00030728  0.99984975]
Iteration 6: [ 1.00009128  2.00002134 -1.00003115  0.9999881 ]
Iteration 7: [ 1.00000836  2.00000117 -1.00000275  0.99999922]
Iteration 8: [ 1.00000067  2.00000002 -1.00000021  0.99999996]
Iteration 9: [ 1.00000004  1.99999999 -1.00000001  1.        ]
Iteration 10: [ 1.  2. -1.  1.]
Solution: [ 1.  2. -1.  1.]
Error: [  2.06480930e-08  -1.25551054e-08   3.61417563e-11   0.00000000e+00]

Программа для решения произвольного номера. уравнений с использованием Matlab

В следующем коде используется формула $$x_{i}^{(k+1)}={\frac {1}{a_{ii}}}\left(b_{i}-\sum _{j<i}a_{ij}x_{j}^{(k+1)}-\sum _{j>i}a_{ij}x_{j}^{(k)}\right),\quad {\begin{array}{l}i=1,2,\ldots ,n\\k=0,1,2,\ldots \end{array}}$$

function x = gauss_seidel(A, b, x, iters)
    for i = 1:iters
        for j = 1:size(A,1)
            x(j) = (b(j) - sum(A(j,:)'.*x) + A(j,j)*x(j)) / A(j,j);
        end
    end
end

См. также

• Метод сопряженных градиентов
• Распространение гауссовских убеждений
• Итерационный метод: Линейные системы
• Метод Качмарца («метод, ориентированный на строки, тогда как метод Гаусса-Зейделя является «ориентированным на столбцы». См., например, эту статью .)
• Расщепление матрицы
• Итерация Ричардсона

Ссылки

• Гаусс, Карл Фридрих (1903), Сочинения (на немецком языке), том. 9, Геттинген: Королевское общество наук .
• Голуб, Джин Х .; Ван Лоан, Чарльз Ф. (1996), Матричные вычисления (3-е изд.), Балтимор: Джонс Хопкинс, ISBN  978-0-8018-5414-9 .
• Блэк, Ноэль и Мур, Ширли. «Метод Гаусса-Зейделя» . Математический мир .

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

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

• «Метод Зейделя» , Математическая энциклопедия , EMS Press , 2001 [1994]
• Гаусса-Зейделя с сайта www.math-linux.com.
• Гаусса – Зейделя из Института холистических численных методов.
• Итерация Гаусса Зиделя с сайта www.geocities.com
• Метод Гаусса-Зейделя
• Биксон
• Код Матлаба
• Пример кода C