Метод Гаусса — Зейделя решения системы линейных уравнений
Метод Гаусса — Зейделя (метод Зейделя, процесс Либмана, метод последовательных замещений) — является классическим итерационным методом решения системы линейных уравнений.
Назван в честь Зейделя и Гаусса.
Постановка задачи
Возьмём систему: [math]\displaystyle{ A\vec{x}=\vec{b} }[/math], где [math]\displaystyle{ 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) }[/math]
Или [math]\displaystyle{ \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. }[/math]
И покажем, как её можно решить с использованием метода Гаусса-Зейделя.
Метод
Чтобы пояснить суть метода, перепишем задачу в виде
- [math]\displaystyle{ \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. }[/math]
Здесь в [math]\displaystyle{ j }[/math]-м уравнении мы перенесли в правую часть все члены, содержащие [math]\displaystyle{ x_i }[/math], для [math]\displaystyle{ i \gt j }[/math]. Эта запись может быть представлена как
- [math]\displaystyle{ (\mathrm{L} + \mathrm{D} )\vec{x} = -\mathrm{U} \, \vec{x} + \vec{b}, }[/math]
где в принятых обозначениях [math]\displaystyle{ \mathrm{D} }[/math] означает матрицу, у которой на главной диагонали стоят соответствующие элементы матрицы [math]\displaystyle{ A }[/math], а все остальные нули; тогда как матрицы [math]\displaystyle{ \mathrm{U} }[/math] и [math]\displaystyle{ \mathrm{L} }[/math] содержат верхнюю и нижнюю треугольные части [math]\displaystyle{ A }[/math], на главной диагонали которых нули.
Итерационный процесс в методе Гаусса — Зейделя строится по формуле
- [math]\displaystyle{ (\mathrm{L} + \mathrm{D})\vec{x}^{(k+1)} = -\mathrm{U} \vec{x}^{(k)} + \vec{b}, \quad k = 0, 1, 2, \ldots }[/math]
после выбора соответствующего начального приближения [math]\displaystyle{ \vec{x}^{(0)} }[/math].
Метод Гаусса — Зейделя можно рассматривать как модификацию метода Якоби. Основная идея модификации состоит в том, что новые значения [math]\displaystyle{ \vec{x}^{(i)} }[/math] используются здесь сразу же по мере получения, в то время как в методе Якоби они не используются до следующей итерации:
- [math]\displaystyle{ \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., }[/math]
где
- [math]\displaystyle{ 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. }[/math]
Таким образом, i-я компонента [math]\displaystyle{ (k+1) }[/math]-го приближения вычисляется по формуле
- [math]\displaystyle{ 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. }[/math]
Например, при [math]\displaystyle{ n=3 }[/math]:
- [math]\displaystyle{ 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 }[/math], то есть [math]\displaystyle{ x_1^{(k+1)}= c_{11}x_{1}^{(k)} + c_{12}x_{2}^{(k)}+ c_{13}x_{3}^{(k)} + d_1, }[/math]
- [math]\displaystyle{ 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 }[/math], то есть [math]\displaystyle{ x_2^{(k+1)}= c_{21}x_{1}^{(k+1)} + c_{22}x_{2}^{(k)}+ c_{23}x_{3}^{(k)} + d_2, }[/math]
- [math]\displaystyle{ 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 }[/math], то есть [math]\displaystyle{ x_3^{(k+1)}= c_{31}x_{1}^{(k+1)} + c_{32}x_{2}^{(k+1)}+ c_{33}x_{3}^{(k)} + d_3. }[/math]
Условие сходимости
Приведём достаточное условие сходимости метода.
Теорема. Пусть [math]\displaystyle{ \| \mathrm{A}_2 \| \lt 1 }[/math], где [math]\displaystyle{ \mathrm{A}_2 = -(\mathrm{L} + \mathrm{D})^{-1} \, \mathrm{U}, \quad (\mathrm{L} + \mathrm{D})^{-1} }[/math] – матрица, обратная к [math]\displaystyle{ (\mathrm{L} + \mathrm{D}) }[/math]. Тогда при любом выборе начального приближения [math]\displaystyle{ \vec{x}^{(0)} }[/math]:
|
Условие окончания
Условие окончания итерационного процесса Зейделя при достижении точности [math]\displaystyle{ \varepsilon }[/math] в упрощённой форме имеет вид:
- [math]\displaystyle{ \parallel x^{(k+1)}-x^{(k)}\parallel \le \varepsilon }[/math]
Более точное условие окончания итерационного процесса имеет вид
- [math]\displaystyle{ \parallel Ax^{(k)}-b\parallel \le \varepsilon }[/math]
и требует больше вычислений. Хорошо подходит для разреженных матриц.
Пример реализации на 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.sqrt(sum((x_new[i] - x[i]) ** 2 for i in range(n))) <= eps
x = x_new
return x