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

Метод Гаусса—Зейделя^[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}=-{\frac {a_{ij}}{a_{ii}}},\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

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

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

Теорема.
Пусть

\|\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)}\!

:

метод Гаусса-Зейделя сходится;
скорость сходимости метода равна скорости сходимости геометрической прогрессии со знаменателем $q=\|\mathrm {A} _{2}\|\!$ ;
верна оценка погрешности: $\|{\vec {x}}^{(k)}-{\vec {x}}\|=q^{k}\,\|{\vec {x}}^{(0)}-{\vec {x}}\|\!$ .

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

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

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

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

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

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

Пример алгоритма на с++

 // Условие сходимости
bool converge(double *xk, double *xkp)
{
    for (int i = 0; i < n; i++) 
    {
        if (fabs(xk[i] - xkp[i]) >= eps) 
            return false;
    }
    return true;
}

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

do
{
    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]);
        p[i] = x[i];
        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;
    s, v, m: double;
begin
  // Проверка на совместность
  for i := 1 to n do
    begin
      s := 0;
      for j := 1 to n do
        if j <> i then s := s + abs(a[i, j]);
      if s >= abs(a[i, i]) then
        begin
          writeln('SLAE is inconsistent');
          exit
        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])-abs(v);
      end;
  until m < e;
  
  // Вывод корней
  writeln('roots: ');
  for i := 1 to n do
    writeln('x[', i, ']= ', x[i]:0:4);
end;

Разбор алгоритма

Суть в том, что стандартная система вида

$~a_{11}c_{1}+a_{12}c_{2}+a_{13}c_{3}=b_{1}$

$~a_{21}c_{1}+a_{22}c_{2}+a_{23}c_{3}=b_{2}$

$~a_{31}c_{1}+a_{32}c_{2}+a_{33}c_{3}=b_{3}$

(где a - вектор коэффициентов, b - вектор свободных членов, а c - вектор неизвестных) преобразуется так, чтобы их могла решать машина.

$~c_{1}={\frac {b_{1}-a_{12}c_{2}-a_{13}c_{3}}{a_{11}}}$

$~c_{2}={\frac {b_{2}-a_{21}c_{1}-a_{23}c_{3}}{a_{22}}}$

$~c_{3}={\frac {b_{3}-a_{31}c_{1}-a_{32}c_{2}}{a_{33}}}$

Перед первой итерацией, вектор (массив) C должен быть инициализирован какими-либо элементами. Тут мнения разделились: одни считают, что предпочтительнее использовать вектор свободных членов (b), другие - нули. От этого выбора зависит только количество итераций (циклов), за которое будет получен ответ с заданной точностью корней (т.е. после каждой итерации Ваши корни будут всё ближе и ближе подбираться к истинным значениям), поэтому для того, чтобы решить, когда пора останавливаться, нужно сравнить текущий вектор С с предыдущим, критерий обычно называют "эпсилон" (и записывают как eps).

#include <math.h>
int main() {
int itr, k, m=3;
float a[m][m], b[m], c[m], c_old[m], temp, eps=0.000001;
// Тут надо вставить код формирования массивов A, B и записать что-нибудь в C
do   // Поехали решать систему ;)
{
    itr++; // Это просто счётчик итераций
    for (k=0; k<m; k++)
        c_old[k]=c[k];
    c[0]=(b[0]-a[0][1]*c[1]-a[0][2]*c[2])/a[0][0];
    c[1]=(b[1]-a[1][0]*c[0]-a[1][2]*c[2])/a[1][1];
    c[2]=(b[2]-a[2][0]*c[0]-a[2][1]*c[1])/a[2][2];
    temp=0;
    for (k=0; k<m; k++)
        temp+=fabs(c_old[k] - c[k]);   // Это мы так вычисляем критерий
} while (temp > eps);
return 0; }

Метод решения СЛАУ был придуман Якоби, но он говорил, что нельзя смешивать старые С и новые, т.е. на i-том шаге вычисляется i-е С по i-1, т.е. по предыдущим. Зейдель с Гауусом показали, что нарушение этого правила не только не приведет к ошибке при нахождении корней, но и может уменьшить количество потребовавшихся итераций для нахождения решения с заданной точностью.

В случае, когда порядок системы большой (много уравнений), человеку становится тяжело выводить каждое уравнение, поэтому он задаёт решение хитрым образом. Если присмотреться к уравнениям, то можно узреть закономерность: для вычисления Cx мы отнимаем от Bx сумму всех (кроме одного) Axy*Cx (т.е. Y пробегает от 0 до m-1, перескакивая через X) и делим всё это на Bxx. В свете последних открытий код выше нужно переписать так:

#include <math.h>
int main() {
int itr, k, l, m=3;
float a[m][m], b[m], c[m], c_old[m], temp, var, eps=0.000001;
// Тут надо вставить код формирования массивов A, B и записать что-нибудь в C
do
{
    itr++;
    temp=0;
    for (k=0; k<m; k++)
    {
        c_old[k]=c[k];
        var=0;
        for (l=0; l<m; l++)
            if (k!=l) var+=a[k][l]*c[l];
        c[k]=(b[k]-var)/a[k][k];
        temp+=fabs(c_old[k]-c[k]);
    }
} while (temp > eps);
return 0; }

Теперь порядок системы (количество уравнений в ней) зависит только от значения переменной m. Компьютер для того и был создан, чтобы делать за людей монотонную работу. :)

Примечания

↑ Людвиг Зейдель (1821—1896) — немецкий астроном и математик, Карл Фридрих Гаусс (1777—1855) — немецкий математик, астроном и физик

См. также

[1] Людвиг Зейдель (1821—1896) — немецкий астроном и математик, Карл Фридрих Гаусс (1777—1855) — немецкий математик, астроном и физик

[1]

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

Содержание

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

Метод

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

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

Пример алгоритма на с++

Пример алгоритма на pascal (для кв. системы)

Разбор алгоритма

Примечания

См. также

Навигация

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

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

Метод

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

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

Пример алгоритма на с++

Пример алгоритма на pascal (для кв. системы)

Разбор алгоритма

Примечания

См. также

Навигация

Поиск