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

(перенаправлено с «Метод Зейделя»)

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

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

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

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

Или  

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

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

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

 

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

 

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

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

 

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

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

 

где

 

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

 

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

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

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

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

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

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

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

 

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

 

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

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

  1 #include <iostream>
  2 #include <cmath>
  3 
  4 using namespace std;
  5 
  6 // Условие окончания
  7 bool converge(double xk[10], double xkp[10], int n, double eps)
  8 {
  9 	double norm = 0;
 10 	for (int i = 0; i < n; i++)
 11 		norm += (xk[i] - xkp[i]) * (xk[i] - xkp[i]);
 12 	return (sqrt(norm) < eps);
 13 }
 14 
 15 double okr(double x, double eps)
 16 {
 17 	int i = 0;
 18 	while (eps != 1)
 19 	{
 20 		i++;
 21 		eps *= 10;
 22 	}
 23 	int okr = pow(double(10), i);
 24 	x = int(x * okr + 0.5) / double(okr);
 25 
 26 	return x;
 27 }
 28 
 29 bool diagonal(double a[10][10], int n)
 30 {
 31 	int i, j, k = 1;
 32 	double sum;
 33 	for (i = 0; i < n; i++) {
 34 		sum = 0;
 35 		for (j = 0; j < n; j++) sum += a[i][j];
 36 		sum -= a[i][i];
 37 		if (sum > a[i][i]) 
 38 		{
 39 			k = 0; 
 40 			cout << a[i][i] << " < " << sum << endl;
 41 		}
 42 		else
 43 		{
 44 			cout << a[i][i] << " > " << sum << endl;
 45 		}
 46 		
 47 
 48 	}
 49 
 50 	return (k == 1);
 51 
 52 }
 53 
 54 
 55 
 56 
 57 int main()
 58 {
 59 	setlocale(LC_ALL, "");
 60 
 61 	double eps, a[10][10], b[10], x[10], p[10];
 62 	int n, i, j, m = 0;
 63 	int method;
 64 	cout << "Введите размер квадратной матрицы: ";
 65 	cin >> n;
 66 	cout << "Введите точность вычислений: ";
 67 	cin >> eps;
 68 	cout << "Заполните матрицу А: " << endl << endl;
 69 	for (i = 0; i < n; i++)
 70 		for (j = 0; j < n; j++)
 71 		{
 72 			cout << "A[" << i << "][" << j << "] = ";
 73 			cin >> a[i][j];
 74 		}
 75 	cout << endl << endl;
 76 	cout << "Ваша матрица А: " << endl << endl;
 77 	for (i = 0; i < n; i++)
 78 	{
 79 		for (j = 0; j < n; j++)
 80 			cout << a[i][j] << " ";
 81 		cout << endl;
 82 	}
 83 
 84 	cout << endl;
 85 
 86 	cout << "Заполните столбец свободных членов: " << endl << endl;
 87 	for (i = 0; i < n; i++)
 88 	{
 89 		cout << "В[" << i + 1 << "] = ";
 90 		cin >> b[i];
 91 	}
 92 
 93 	cout << endl << endl;
 94 
 95 
 96 	/*
 97 	Ход метода, где:
 98 	a[n][n] - Матрица коэффициентов
 99 	x[n], p[n] - Текущее и предыдущее решения
100 	b[n] - Столбец правых частей
101 	Все перечисленные массивы вещественные и
102 	должны быть определены в основной программе,
103 	также в массив x[n] следует поместить начальное
104 	приближение столбца решений (например, все нули)
105 	*/
106 
107 	for (int i = 0; i < n; i++)
108 		x[i] = 1;
109 
110 	cout << "Диагональное преобладание: " << endl;
111 	if (diagonal(a, n)) {
112 		do
113 		{
114 			for (int i = 0; i < n; i++)
115 				p[i] = x[i];
116 			for (int i = 0; i < n; i++)
117 			{
118 				double var = 0;
119 				for (int j = 0; j < i; j++)
120 					var += (a[i][j] * x[j]);
121 				for (int j = i + 1; j < n; j++)
122 					var += (a[i][j] * p[j]);
123 				x[i] = (b[i] - var) / a[i][i];
124 			}
125 			m++;
126 		} while (!converge(x, p, n, eps));
127 
128 
129 
130 		cout << "Решение системы:" << endl << endl;
131 		for (i = 0; i < n; i++) cout << "x" << i << " = " << okr(x[i], eps) << "" << endl;
132 		cout << "Итераций: " << m << endl;
133 	}
134 	else {
135 		cout << "Не выполняется преобладание диагоналей" << endl;
136 	}
137 
138 
139 	system("pause");
140 	return 0;
141 }

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

from math import sqrt
import numpy as np

def seidel(A, b, eps):
    n = len(A)
    x = [.0 for i in range(n)]

    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 = sqrt(sum((x_new[i] - x[i]) ** 2 for i in range(n))) <= eps
        x = x_new

    return x

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

ПримечанияПравить

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