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

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

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

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

Возьмём систему: , где

Или

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

Метод

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

Здесь в -м уравнении мы перенесли в правую часть все члены, содержащие , для . Эта запись может быть представлена как

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

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

после выбора соответствующего начального приближения .

Метод Гаусса — Зейделя можно рассматривать как модификацию метода Якоби. Основная идея модификации состоит в том, что новые значения используются здесь сразу же по мере получения, в то время как в методе Якоби они не используются до следующей итерации:

где

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

Например, при :

, то есть
, то есть
, то есть

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

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

Теорема.
Пусть , где – матрица, обратная к . Тогда при любом выборе начального приближения :
  1. метод Гаусса-Зейделя сходится;
  2. скорость сходимости метода равна скорости сходимости геометрической прогрессии со знаменателем ;
  3. верна оценка погрешности: .

Условие окончания

Условие окончания итерационного процесса Зейделя при достижении точности в упрощённой форме имеет вид:

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

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

Пример реализации на 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

См. также

Примечания

Ссылки