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

Материал из Википедии — свободной энциклопедии
Перейти к: навигация, поиск

Метод Гаусса—Зейделя[1] является классическим итерационным методом решения системы линейных уравнений.

Постановка задачи[править | править вики-текст]

Возьмём систему: 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-м уравнении мы перенесли в правую часть все члены, содержащие x_i , для i > j. Эта запись может быть представлена:

(\mathrm{L} + \mathrm{D} )\vec{x} = -\mathrm{U} \, \vec{x} + \vec{b},

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

Итерационный процесс в методе Гаусса-Зейделя строится по формуле (\mathrm{L} + \mathrm{D} )\vec{x}^{(k+1)} = -\mathrm{U} \vec{x}^{(k)} + \vec{b} ,\quad k = 0, 1, 2, \ldots после выбора соответствующего начального приближения \vec{x}^{(0)}.

Метод Гаусса-Зейделя можно рассматривать как модификацию метода Якоби. Основная идея модификации состоит в том, что новые значения \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)-го приближения вычисляется по формуле:

x_i^{(k+1)}=\sum_{j=1}^{i-1}c_{ij}x_{j}^{(k+1)}+\sum_{j=i+1}^{n}c_{ij}x_{j}^{(k)}+d_i, \quad i=1,\ldots,n

Условие сходимости[править | править вики-текст]

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

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

Условие окончания[править | править вики-текст]

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

\parallel x^{(k+1)}-x^{(k)}\parallel \le \varepsilon

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

\parallel Ax^{(k)}-b\parallel \le \varepsilon

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

Пример реализации на C++[править | править вики-текст]

// Условие окончания
bool converge(double *xk, double *xkp)
{
    double norm = 0;
    for (int i = 0; i < n; i++) 
    {
      norm += (xk[i] - xkp[i])*(xk[i] - xkp[i]);
    }
    if(sqrt(norm) >= eps)
      return false;
    return true;
}
 
/*
    Ход метода, где:
    a[n][n] - Матрица коэффициентов
    x[n], p[n] - Текущее и предыдущее решения
    b[n] - Столбец правых частей
    Все перечисленные массивы вещественные и
    должны быть определены в основной программе,
    также в массив x[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 < i; j++)
            var += (a[i][j] * x[j]);
        for (int j = i; j < n; j++)
            var += (a[i][j] * p[j]);
        x[i] = (b[i] - var) / a[i][i];
    }
}
while (!converge(x, p));

Пример реализации на Pascal[править | править вики-текст]

type ar2d = array [1..50, 1..50] of double;
     ar1d = array [1..50] of double;
 
procedure seidel(n: byte; e: extended; a: ar2d; b: ar1d; x: ar1d);
var i, j: longint;
    v, m: double;
begin
  // Проверка на совместность
  for i := 1 to n do
    begin
      for j := 1 to n do
        if j <> i then
        begin
          if abs(a[i, j]) >= abs(a[i, i]) then
          begin
            writeln('SLAE is inconsistent');
            exit;
          end;
        end;
 
    end;
 
  // Сам алгоритм
  repeat
    m := 0;
    for i := 1 to n do
      begin
        s := 0;
        for j := 1 to n do
          if i <> j then s := s + a[i, j] * x[j];
        v := x[i];
        x[i] := (b[i] - s) / a[i, i];
        m:=abs(x[i]-v);
      end;
  until m < e;
 
  // Вывод корней
  writeln('roots: ');
  for i := 1 to n do
    writeln('x[', i, ']= ', x[i]:0:4);
end;


Пример реализации на Python[править | править вики-текст]

def seidel(m, b, eps):
    n = len(m)
    r = range(n)
    x = [0] * len(m)
    conv = False
    while not conv:
        p = x.copy()
        for i in r:
            var = sum(m[i][j] * x[j] for j in range(i))
            var += sum(m[i][j] * p[j] for j in r)
            x[i] = (b[i] - var) / m[i][i]
 
        conv = sum((x[i]-p[i])**2 for i in r) < eps
    return x

Примечания[править | править вики-текст]

См. также[править | править вики-текст]