Содержание

Решение систем линейных алгебраических уравнений

Данная статья посвящена нахождению корней систем линейных алгебраических уравнений. Методы численного решения таких систем подразделяются на два типа: прямые (конечные) и итерационные (бесконечные). Оба типа полезны и удобны для практических вычислений и каждый из них имеет свои преимущества и недостатки.

Метод исключений (метод Гаусса)

Данный метод является наиболее известным и широко применяемым. Рассмотрим систему из n уравнений с n неизвестными. Обозначим неизвестные через x[1], x[2], … , x[n] и запишем систему в следующем виде:

     a[1,1]*x[1] + a[1,2]*x[2] + ... + a[1,n]*x[n] = b[1]
     a[2,1]*x[1] + a[2,2]*x[2] + ... + a[2,n]*x[n] = b[2]
     ....................................................
     a[n,1]*x[1] + a[n,2]*x[2] + ... + a[n,n]*x[n] = b[n]

Предполагается, что в силу расположения уравнений, a[i,i] <> 0. Если это не так, то, меняя уравнения местами, добиваемся выполнения этого условия. Введем n-1 множителей:

     m[i] = a[i,1]/a[1,1]   , где i = 2, 3, ... , n

и вычтем из каждого i-го уравнения первое, умноженное на m, обозначая

     a'[i,j] = a[i,j] - m[i]*a[1,j]
     b'[i] = - m[i]*b[i]
     i = 2, 3, ... , n;  j = 1, 2, ... , n.

Для всех уравнений, начиная со второго, получим:

     a[i,1] = 0   , где i = 2, 3, ... , n.

Преобразованная система уравнений запишется в следующем виде:

     a[1,1]*x[1] + a[1,2]*x[2] + ... + a[1,n]*x[n] = b[1]
        0   +     a'[2,2]*x[2] + ... + a'[2,n]*x[n] = b'[2]
     ....................................................
        0   +     a'[n,2]*x[2] + ... + a'[n,n]*x[n] = b'[n]

Продолжая таким же образом, исключаем x[2] из последних n-2 уравнений, затем x[3] из последних n-3 и т.д. На некотором k-м этапе мы исключаем x[k] с помощью множителей:

     m^(k-1)_[i] = a^(k-1)_[i,k]/a^(k-1)_[k,k]   , i = k+1, ... , n,

, незабывая, что a(k-1)[k,k] <> 0.

Примечание

^(k-1) обозначает на k-1-й итерации, _[i,k] - индекс эл-та. Тогда

     a^(k)_[i,j] = a^(k-1)_[i,j] - m^(k-1)_[i]*a^(k-1)_[k,j]
     b^(k)_[i] = - m^(k-1)_[i]*b^(k-1)_[k]
     i = k+1, k+2, ... , n;  j = k, k+1, ... , n.

Окончательно, треугольная система уравнений записывается:

     a[1,1]*x[1] + a[1,2]*x[2] + ... + a[1,n]*x[n] = b[1]
                   a[2,2]*x[2] + ... + a[2,n]*x[n] = b[2]
     ....................................................
                                       a[n,n]*x[n] = b[n]

Обратная подстановка для нахождения значений неизвестных задается формулами:

     x[n] = b^(n-1)_[n] / a^(n-1)_[n,n]
     x[n-1] = (b^(n-2)_[n-1] - a^(n-2)_[n-1,n]*x[n]) / a^(n-2)_[n-1,n-1]
     ...................................................................
     x[j] = ( b^(j-1)_[j] - a^(j-1)_[j,n]*x[n] - ... -
            - a^(j-1)_[j,j+1]*x[j+1] ) / a^(j-1)_[j,j]

     j = n-2, n-3, ... , 1.

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

     | a^(k-1)_[k,k] | >= | a^(k-1)_[i,k] |

Поэтому перед началом процесса исключения очередного неизвестного необходимо переставить уравнения так, чтобы максимальный по модулю коэффициент при x[k] попал на главную диагональ (данный способ решения часто называется методом главного элемента).

Итерационный метод Гаусса-Зейделя

Данный метод отличается простотой и легкостью программирования, малой ошибкой округления, но сходится при выполнении некоторых условий.
Рассмотрим все ту же систему из n уравнений с n неизвестными.
По прежнему полагаем, что a[i,i] <>0 для всех i. Тогда k-е приближение к решению будет задаваться формулой:

  x^(k)_[i] = ( b[i] - a[i,1]*x^(k)_[1] - ... - a[i,j-1]*x^(k)_[i-1] - a[i,j+1]*x^(k-1)_[i+1] - a[i,n]*x^(k-1)_[n] ) / a[i,i]
  i = 1, 2, ... , n.

Как известно, любой итерационный процесс, как и настоящий мужик, должен имееть свой конец :) Поэтому вычисления нужно продолжать до тех пор, все x^(k)_[i] не станут достаточно близки к x(k-1)[i], т.е. пока для всех i не будет выполнятся неравенство

     max| x^(k)_[i] - x^(k-1)_[i] | < Eps

, где Eps - некоторое положительное число, задающее точность вычислений. Итерационный метод Гаусса-Зейделя сходится, если удовлетворяется достаточное условие сходимости: для всех i модуль диагонального коэффициента должен быть не меньше суммы модулей остальных коэффициентов данной строки

     | a[i,i] | >= | a[i,1] | + | a[i,2] | + ... + | a[i,j-1] | + | a[i,j+1] | + ... + | a[i,n] | ,

а хотябы для обной строки это неравенство должно выполнятся строго

     | a[i,i] | > | a[i,1] | + | a[i,2] | + ... + | a[i,j-1] | + | a[i,j+1] | + ... + | a[i,n] |.

Сравнение методов

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

метода исключений время пропорционально n3. Если для решения системы требуется меньше чем n итераций, то общие затраты машинного времени будут меньше.

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

А вот и исходник по теме:

{$N+,G+,X+}
 
const
  Comments: boolean = false; { нужно ли печатать пошаговые рез-ты расчета }
  Eps = 0.00001;        { Точность вычислений }
  n_max = 4;            { Кол-во ур-й и неизвестных }
  const_AnB: array[1..n_max,1..n_max+1] of double =
 (( -8,   2,  17,  -5, -119.97),
 { -8*X[1] + 2*X[2] + 17*X[3] - 5*X[4] = -119.97  и т.д. ... }
  (  4, -22,   6,   5,   55.73),
  ( 15,   3,  -5,  -5,   18.77),
  ( -4,  -4,   5,  14,  -79.42));
{ Описание типов матрицы уравнеиний и массива найденных Х }
Type TMatr = array[1..n_max,1..n_max+1] of double;
Type TXMarr = array[1..n_max] of double;
 
procedure PrintX(const X: TXMarr; n: byte; cnt: word);
{ Печатает найденные Х }
var
  k,       { Позиция заданного корня в массиве }
  i: byte; { Номер текущего корня для вывода }
begin
  if n = 0 then
    exit;
  k:=low(X);
  if cnt = 0 then
    begin
      write(' N       x1');
      for i:=2 to n do
        write('         x',i);
      writeln;
    end;
  write(cnt:2);
  for i:=1 to n do
    begin
      Write('   ',X[k]:3:5);
      inc(k);
    end;
  writeln;
{
  При решении систем, с большим кол-вом неизвестных, возможно,
  прийдется переделать эту процедуру для корректного вывода
  всех решений.
}
end; { PrintX }
 
procedure PrintMatr(var AnB: TMatr; n: byte);
{ Печатает заданную матрицу }
var
  i,j: byte;
begin
  for i:=1 to n do
    begin
      for j:=1 to n+1 do
        write(AnB[i,j]:3:5,'   ');
      writeln;
    end;
  writeln;
{
  При решении систем, с большим кол-вом неизвестных, возможно,
  прийдется переделать эту процедуру для корректного вывода
  всех решений.
}
end; { PrintMatr }
 
function Normalize(var AnB: TMatr; n: byte): byte;
{
 Располагает на главной диагонали наибольшие элементы в столбцах
 Если один из них равен нулю,то систему решить не получится ...
 Возвращает 0 если систему заданными методами решить не получится,
 1 - с-му можно решить только методом Гаусса
 2 - с-му можно решать любым методом
}
var
  max: double; { Переменная для поиска макс. эл-та }
  imax: byte;  { Переменная для поиска макс. эл-та }
  tmp: double; { Переменная для обмена элементов }
  i,j: byte;   { Счетчики циклов }
begin
  Normalize:=2;
  for j:=1 to n do
    begin
      max:=Abs(AnB[1,j]);
      imax:=1;
      for i:=2 to n do
        if Abs(AnB[i,j]) > max then
          begin
            max:=Abs(AnB[i,j]);
            imax:=i;
          end;
      if max < Eps then
        begin
          Normalize:=0;
          Writeln('Error: a[i,i]=0!');
          exit;
        end
      else
        if j <> imax then
          for i:=1 to n+1 do
            begin
              tmp:=AnB[j,i];
              AnB[j,i]:=AnB[imax,i];
              AnB[imax,i]:=tmp;
            end;
    end;
  for i:=1 to n do
    begin
      tmp:=0;
      for j:=1 to n do
        if i <> j then
          tmp:=tmp+Abs(AnB[i,j]);
      if Abs(AnB[i,i]) < tmp then
        Normalize:=1;
    end;
end; { Normalize }
 
procedure Metod1(const Matr; n: byte);
var
  k    : byte;
  M_   : TMatr absolute Matr;
  AnB  : TMatr;
  X,M  : TXMarr;
function Calc(k: byte): boolean;
var
  i,j  : byte;
begin
  Calc:= true;
  for i:=k+1 to n do
    begin
      M[i]:= AnB[i,k]/AnB[k,k];
      if M[i] > 1 then
        begin
          Calc:= false;
          exit;
        end;
    end;
 
  for i:=k+1 to n do
    for j:=k to n+1 do
      AnB[i,j]:=AnB[i,j]-M[i]*AnB[k,j];
 
  if Comments then
    PrintMatr(AnB,n);
end; { Calc }
function Summ(j: byte): double;
var
  i : byte;
  res : double;
begin
  res:=AnB[j,n+1];
  for i:=n downto j+1 do
    res:=res - AnB[j,i]*X[i];
  Summ:=res;
end; { Summ }
begin
  Writeln('Метод Гаусса:');
  AnB:=M_;
  if Normalize(AnB,n) = 0 then
    begin
      writeln('Error: Невозможно решить систему уравнений!');
      exit;
    end;
 
  for k:=1 to n-1 do
    if not Calc(k) then
      begin
        Writeln('Error: Невыполняется условие | a^(k-1)_[k,k] | >= | a^(k-1)_[i,k] | !');
        exit;
      end;
 
  X[n]:=AnB[n,n+1]/AnB[n,n];
  for k:=n-1 downto 1 do
    X[k]:=Summ(k)/AnB[k,k];
 
  PrintX(X,n,0);
  Writeln;
end; { Metod1 }
 
procedure Metod2(const Matr; n: byte);
var
  M_   : TMatr absolute Matr;
  AnB  : TMatr;
  X    : TXMarr;
  i,j: byte;
  tmp: double;
  delta: double;
  cnt: word;
function Summ(i: byte): double;
var
  j: byte;
  res: double;
begin
  res:=AnB[i,n+1];
  for j:=1 to n do
    if i <> j then
      res:=res-AnB[i,j]*X[j];
  Summ:=res;
end; { Summ }
begin
  writeln('Метод Гаусса-Зейделя:');
  AnB:=M_;
  For i:=1 to n do
    X[i]:=0;
  if Normalize(AnB,n) <> 2 then
    begin
      writeln('Error: !');
      exit;
    end;
  delta:=1;
  cnt:=0;
  while delta > Eps do
    begin
      if Comments then
        PrintX(X,n,cnt);
      delta:=0;
      for i:=1 to n do
        begin
          tmp:=Summ(i)/AnB[i,i];
          if delta < Abs(tmp-X[i]) then
            delta:=Abs(tmp-X[i]);
          X[i]:=tmp;
        end;
      inc(cnt);
    end;
 
  PrintX(X,n,cnt);
 
  writeln('Заданная точность вычислений достигнута на ',cnt,' шаге.');
  Writeln;
end; { Metod2 }
 
begin
  Writeln(#10,'Численное решение систем линейных алгебраических уравнений:',#10);
  Metod1(const_AnB,n_max);
  Writeln('Press Enter to continue');
  readln;
  Metod2(const_AnB,n_max);
  Writeln('Press Enter if you wanna DOS ;)');
  readln;
end.

Автор: Labinskiy Nikolay aka e-moe © 2005