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

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

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

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

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

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

Или

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

Метод[править | править код]

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

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

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

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

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

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

где

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

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

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

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

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

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

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

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

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

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

Пример реализации на 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 + 1; 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[править | править код]

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

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

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