Метод Гаусса — Зейделя решения системы линейных уравнений

Не следует путать с методом покоординатного спуска Гаусса — Зейделя — методом поиска экстремума.

Метод Гаусса — Зейделя (метод Зейделя, процесс Либмана, метод последовательных замещений) — является классическим итерационным методом решения системы линейных уравнений.

Назван в честь Зейделя и Гаусса.

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

Возьмём систему: $A\stackrel{\to <!--\; \to \; -->}{x}=\stackrel{\to <!--\; \to \; -->}{b}$$A\vec{x}=\vec{b}$ , где $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)$ 

Или $\left\{ \begin{array}{rcl} a_{11}x_1 + \ldots + a_{1n}x_n& = & b_{1} \\ & &\\ a_{n1}x_1 + \ldots + a_{nn}x_n & = & b_{n} \end{array} \right.$ 

И покажем, как её можно решить с использованием метода Гаусса-Зейделя.

МетодПравить

Чтобы пояснить суть метода, перепишем задачу в виде

$\left\{ \begin{array}{lcr} a_{11}x_1 & = &-a_{12}x_2 - a_{13}x_3 -\ldots - a_{1n}x_n + b_1\\ a_{21}x_1 + a_{22}x_2 & = & -a_{23}x_3 - \ldots - a_{2n}x_n + b_2\\ \ldots & &\\ a_{(n-1)1}x_1 + a_{(n-1)2}x_2 +\ldots+ a_{(n-1)(n-1)}x_{n-1} & = & -a_{(n-1)n}x_n + b_{n-1}\\ a_{n1}x_1 + a_{n2}x_2 +\ldots+ a_{n(n-1)}x_{n-1}+ a_{nn}x_n & = & b_n \end{array} \right.$ 

Здесь в $j$$j$ -м уравнении мы перенесли в правую часть все члены, содержащие $x_i$$x_{i}$ , для $i>j$$i>j$ . Эта запись может быть представлена как

$(\mathrm{L}+\mathrm{D})\stackrel{\to <!--\; \to \; -->}{x}=-<!--\; -\; -->\mathrm{U}\stackrel{\to <!--\; \to \; -->}{x}+\stackrel{\to <!--\; \to \; -->}{b},$$(\mathrm{L} + \mathrm{D} )\vec{x} = -\mathrm{U} \, \vec{x} + \vec{b},$ 

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

Итерационный процесс в методе Гаусса — Зейделя строится по формуле

$(\mathrm{L}+\mathrm{D}){\stackrel{\to <!--\; \to \; -->}{x}}^{(k+1)}=-<!--\; -\; -->\mathrm{U}{\stackrel{\to <!--\; \to \; -->}{x}}^{(k)}+\stackrel{\to <!--\; \to \; -->}{b},k=0,1,2,\dots <!--\; \dots \; -->$$(\mathrm {L} +\mathrm {D} ){\vec {x}}^{(k+1)}=-\mathrm {U} {\vec {x}}^{(k)}+{\vec {b}},\quad k=0,1,2,\ldots$ 

после выбора соответствующего начального приближения ${\stackrel{\to <!--\; \to \; -->}{x}}^{(0)}$$\vec{x}^{(0)}$ .

Метод Гаусса — Зейделя можно рассматривать как модификацию метода Якоби. Основная идея модификации состоит в том, что новые значения ${\stackrel{\to <!--\; \to \; -->}{x}}^{(i)}$$\vec{x}^{(i)}$  используются здесь сразу же по мере получения, в то время как в методе Якоби они не используются до следующей итерации:

$\left\{\begin{array}{ccccccccccc} {x_{1}}^{(k+1)} &=& c_{12}{x_2^{(k)}} &+& c_{13}x_3^{(k)}&+& {\ldots}&+& c_{1n}{x_n}^{(k)} &+& d_1 \\ {x_{2}}^{(k+1)} &=& c_{21}{x_1^{(k+1)}} &+& c_{23}x_3^{(k)}&+& {\ldots}&+& c_{2n}{x_n}^{(k)} &+& d_2 \\ \ldots & & & & & & & & & & \\ {x_{n}}^{(k+1)} &=& c_{n1}{x_1^{(k+1)}} &+& c_{n2}{x_2^{(k+1)}}&+& {\ldots}&+& c_{n(n-1)}{x_{n-1}}^{(k+1)} &+& d_n \end{array}\right.,$ 

где

$c_{ij}={\begin{cases}-{\frac {a_{ij}}{a_{ii}}},&j\neq i,\\0,&j=i.\end{cases}}\quad d_{i}={\frac {b_{i}}{a_{ii}}},\quad i=1,\ldots ,n.$ 

Таким образом, i-я компонента $(k+1)$$(k+1)$ -го приближения вычисляется по формуле

$x_i^{(k+1)}=\underset{j=1}{\overset{i-<!--\; -\; -->1}{\sum <!--\; \sum \; -->}}{c}_{ij}{x}_j^{(k+1)}+\underset{j=i}{\overset{n}{\sum <!--\; \sum \; -->}}{c}_{ij}{x}_j^{(k)}+d_i,i=1,\dots <!--\; \dots \; -->,n.$$x_{i}^{(k+1)}=\sum _{j=1}^{i-1}c_{ij}x_{j}^{(k+1)}+\sum _{j=i}^{n}c_{ij}x_{j}^{(k)}+d_{i},\quad i=1,\ldots ,n.$ 

Например, при $n=3$$n=3$ :

$x_1^{(k+1)}=\underset{j=1}{\overset{1-<!--\; -\; -->1}{\sum <!--\; \sum \; -->}}{c}_{1j}{x}_j^{(k+1)}+\underset{j=1}{\overset{3}{\sum <!--\; \sum \; -->}}{c}_{1j}{x}_j^{(k)}+d_1$$x_{1}^{{(k+1)}}=\sum _{{j=1}}^{{1-1}}c_{{1j}}x_{{j}}^{{(k+1)}}+\sum _{{j=1}}^{{3}}c_{{1j}}x_{{j}}^{{(k)}}+d_{1}$ , то есть $x_1^{(k+1)}=c_{11}{x}_1^{(k)}+c_{12}{x}_2^{(k)}+c_{13}{x}_3^{(k)}+d_1,$$x_{1}^{(k+1)}=c_{11}x_{1}^{(k)}+c_{12}x_{2}^{(k)}+c_{13}x_{3}^{(k)}+d_{1},$ 

$x_2^{(k+1)}=\underset{j=1}{\overset{2-<!--\; -\; -->1}{\sum <!--\; \sum \; -->}}{c}_{2j}{x}_j^{(k+1)}+\underset{j=2}{\overset{3}{\sum <!--\; \sum \; -->}}{c}_{2j}{x}_j^{(k)}+d_2$$x_{2}^{{(k+1)}}=\sum _{{j=1}}^{{2-1}}c_{{2j}}x_{{j}}^{{(k+1)}}+\sum _{{j=2}}^{{3}}c_{{2j}}x_{{j}}^{{(k)}}+d_{2}$ , то есть $x_2^{(k+1)}=c_{21}{x}_1^{(k+1)}+c_{22}{x}_2^{(k)}+c_{23}{x}_3^{(k)}+d_2,$$x_{2}^{(k+1)}=c_{21}x_{1}^{(k+1)}+c_{22}x_{2}^{(k)}+c_{23}x_{3}^{(k)}+d_{2},$ 

$x_3^{(k+1)}=\underset{j=1}{\overset{3-<!--\; -\; -->1}{\sum <!--\; \sum \; -->}}{c}_{3j}{x}_j^{(k+1)}+\underset{j=3}{\overset{3}{\sum <!--\; \sum \; -->}}{c}_{3j}{x}_j^{(k)}+d_3$$x_{3}^{{(k+1)}}=\sum _{{j=1}}^{{3-1}}c_{{3j}}x_{{j}}^{{(k+1)}}+\sum _{{j=3}}^{{3}}c_{{3j}}x_{{j}}^{{(k)}}+d_{3}$ , то есть $x_3^{(k+1)}=c_{31}{x}_1^{(k+1)}+c_{32}{x}_2^{(k+1)}+c_{33}{x}_3^{(k)}+d_3.$$x_{3}^{(k+1)}=c_{31}x_{1}^{(k+1)}+c_{32}x_{2}^{(k+1)}+c_{33}x_{3}^{(k)}+d_{3}.$ 

Условие сходимостиПравить

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

Теорема. Пусть $\|\mathrm {A} _{2}\|<1$ , где ${\mathrm{A}}_2=-<!--\; -\; -->(\mathrm{L}+\mathrm{D}{)}^{-<!--\; -\; -->1}\mathrm{U},(\mathrm{L}+\mathrm{D}{)}^{-<!--\; -\; -->1}$$\mathrm {A} _{2}=-(\mathrm {L} +\mathrm {D} )^{-1}\,\mathrm {U} ,\quad (\mathrm {L} +\mathrm {D} )^{-1}$  – матрица, обратная к $(\mathrm{L}+\mathrm{D})$$(\mathrm {L} +\mathrm {D} )$ . Тогда при любом выборе начального приближения ${\stackrel{\to <!--\; \to \; -->}{x}}^{(0)}$$\vec{x}^{(0)}$ : метод Гаусса-Зейделя сходится; скорость сходимости метода равна скорости сходимости геометрической прогрессии со знаменателем $q=\Vert <!--\; \Vert \; -->{\mathrm{A}}_2\Vert <!--\; \Vert \; -->$$q=\|\mathrm {A} _{2}\|$ ; верна оценка погрешности: $\Vert <!--\; \Vert \; -->{\stackrel{\to <!--\; \to \; -->}{x}}^{(k)}-<!--\; -\; -->\stackrel{\to <!--\; \to \; -->}{x}\Vert <!--\; \Vert \; -->=q^k\Vert <!--\; \Vert \; -->{\stackrel{\to <!--\; \to \; -->}{x}}^{(0)}-<!--\; -\; -->\stackrel{\to <!--\; \to \; -->}{x}\Vert <!--\; \Vert \; -->$$\|{\vec {x}}^{(k)}-{\vec {x}}\|=q^{k}\,\|{\vec {x}}^{(0)}-{\vec {x}}\|$ .

Условие окончанияПравить

Условие окончания итерационного процесса Зейделя при достижении точности $\mathrm{\epsilon <!--\; \epsilon \; -->}$$\varepsilon$  в упрощённой форме имеет вид:

$\parallel <!--\; \parallel \; -->x^{(k+1)}-<!--\; -\; -->x^{(k)}\parallel <!--\; \parallel \; -->\le <!--\; \le \; -->\mathrm{\epsilon <!--\; \epsilon \; -->}$$\parallel x^{(k+1)}-x^{(k)}\parallel \le \varepsilon$ 

Более точное условие окончания итерационного процесса имеет вид

$\parallel <!--\; \parallel \; -->A{x}^{(k)}-<!--\; -\; -->b\parallel <!--\; \parallel \; -->\le <!--\; \le \; -->\mathrm{\epsilon <!--\; \epsilon \; -->}$$\parallel Ax^{(k)}-b\parallel \le \varepsilon$ 

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

Пример реализации на C++Править

#include <iostream>
#include <cmath>

using namespace std;

// Условие окончания
bool converge(double xk[10], double xkp[10], int n, double eps)
{
	double norm = 0;
	for (int i = 0; i < n; i++)
		norm += (xk[i] - xkp[i]) * (xk[i] - xkp[i]);
	return (sqrt(norm) < eps);
}

double okr(double x, double eps)
{
	int i = 0;
	double neweps = eps;
	while (neweps < 1)
	{
		i++;
		neweps *= 10;
	}
	int okr = pow(double(10), i);
	x = int(x * okr + 0.5) / double(okr);

	return x;
}

bool diagonal(double a[10][10], int n)
{
	int i, j, k = 1;
	double sum;
	for (i = 0; i < n; i++) {
		sum = 0;
		for (j = 0; j < n; j++) sum += abs(a[i][j]);
		sum -= abs(a[i][i]);
		if (sum > a[i][i]) 
		{
			k = 0; 
			cout << a[i][i] << " < " << sum << endl;
		}
		else
		{
			cout << a[i][i] << " > " << sum << endl;
		}
		

	}

	return (k == 1);

}




int main()
{
	setlocale(LC_ALL, "");

	double eps, a[10][10], b[10], x[10], p[10];
	int n, i, j, m = 0;
	int method;
	cout << "Введите размер квадратной матрицы: ";
	cin >> n;
	cout << "Введите точность вычислений: ";
	cin >> eps;
	cout << "Заполните матрицу А: " << endl << endl;
	for (i = 0; i < n; i++)
		for (j = 0; j < n; j++)
		{
			cout << "A[" << i << "][" << j << "] = ";
			cin >> a[i][j];
		}
	cout << endl << endl;
	cout << "Ваша матрица А: " << endl << endl;
	for (i = 0; i < n; i++)
	{
		for (j = 0; j < n; j++)
			cout << a[i][j] << " ";
		cout << endl;
	}

	cout << endl;

	cout << "Заполните столбец свободных членов: " << endl << endl;
	for (i = 0; i < n; i++)
	{
		cout << "В[" << i + 1 << "] = ";
		cin >> b[i];
	}

	cout << endl << endl;


	/*
	Ход метода, где:
	a[n][n] - Матрица коэффициентов
	x[n], p[n] - Текущее и предыдущее решения
	b[n] - Столбец правых частей
	Все перечисленные массивы вещественные и
	должны быть определены в основной программе,
	также в массив x[n] следует поместить начальное
	приближение столбца решений (например, все нули)
	*/

	for (int i = 0; i < n; i++)
		x[i] = 1;

	cout << "Диагональное преобладание: " << endl;
	if (diagonal(a, n)) {
		do
		{
			for (int i = 0; i < n; i++)
				p[i] = x[i];
			for (int i = 0; i < n; i++)
			{
				double var = 0;
				for (int j = 0; j < n; j++)
                    if(j!=i) var += (a[i][j] * x[j]);
				
				x[i] = (b[i] - var) / a[i][i];
			}
			m++;
		} while (!converge(x, p, n, eps));



		cout << "Решение системы:" << endl << endl;
		for (i = 0; i < n; i++) cout << "x" << i << " = " << okr(x[i], eps) << "" << endl;
		cout << "Итераций: " << m << endl;
	}
	else {
		cout << "Не выполняется преобладание диагоналей" << endl;
	}

	system("pause");
	return 0;
}

Пример реализации на PythonПравить

import numpy as np

def seidel(A, b, eps):
    n = len(A)
    x = np.zeros(n)  # zero vector

    converge = False
    while not converge:
        x_new = np.copy(x)
        for i in range(n):
            s1 = sum(A[i][j] * x_new[j] for j in range(i))
            s2 = sum(A[i][j] * x[j] for j in range(i + 1, n))
            x_new[i] = (b[i] - s1 - s2) / A[i][i]

        converge = np.linalg.norm(x_new - x) <= eps
        x = x_new

    return x

См. такжеПравить

• Геометрическая прогрессия
• Метод простой итерации
• Метод Якоби
• Теорема Банаха
• Метод итерации

СсылкиПравить

• https://www.maa.org/press/periodicals/loci/joma/iterative-methods-for-solving-iaxi-ibi-gauss-seidel-method
• https://www3.nd.edu/~zxu2/acms40390F12/Lec-7.3.pdf
• http://mathworld.wolfram.com/Gauss-SeidelMethod.html