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

      x₃ = b₃₁x₁ + b₃₂x₂ +  … + b_3,n–1x_n–1 + b_3nx_n+ c₃ ,

      .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .

      x_n = b_n1x₁ + b_n2x₂ + b_n3x₃ + … + b_n,n–1x_n–1 +  c_n ,

в которой  на главной диагонали матрицы  B находятся нулевые элементы. Остальные элементы выражаются по формулам

      b_ij = –a_ij / a_ii, c_i = b_i / a_ii (i, j = 1, 2, …, n, j ≠ i)

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

1.5.2. Описание метода

     Введем  нижнюю и верхнюю треугольные  матрицы

  0 0 0 … 0  0 b₁₂ b₁₃ … b_1n

            b₂₁ 0 0 … 0  0 0 b₂₃ … b_2n

      B₁ = b₃₁ b₃₂ 0 … 0 ,  B₂ =  0 0 0 … b_3n

            .   .   .   .   .   .   .  .   .   .   .   .   .   .

  b_n1 b_n2 b_n3 … 0  0 0 0 … 0

     Заметим, что B = B₁ + B₂ и поэтому решение x исходной системы удовлетворяет равенству

      x = B₁x + B₂ x + c .

     Выберем начальное приближение x⁽⁰⁾ = [x₁⁽⁰⁾, x₂⁽⁰⁾, …, x_n⁽⁰⁾]^T. Подставляя его в правую часть равенства при верхней треугольной матрице B₂ и вычисляя полученное выражение, находим первое приближение

      x⁽¹⁾ = B₁x⁽⁰⁾ + B₂x⁽¹⁾

Подставляя  приближение x⁽¹⁾, получим

      x⁽²⁾ = B₁x⁽¹⁾ + B₂x⁽²⁾

Продолжая этот процесс далее, получим последовательность x⁽⁰⁾, x⁽¹⁾, …, x⁽ⁿ⁾, приближений к вычисляемых по формуле

      x^(k+1) = B₁^(k+1) + B₂^(k) + c

или в  развернутой форме записи

      x₁^(k+1) =  b₁₂x₂^(k) + b₁₃x₂^(k) + … + b_1nx_n^(k) + c₁ ,

      x₂^(k+1) = b₂₁x₁^(k+1) +  b₂₃x₃^(k) +  … + b_2nx_n^(k) + c₂ ,

      x₃^(k+1) = b₃₁x₁^(k+1) + b₃₂x₂^(k+1) +  … + b_3nx_n^(k) + c₃ ,

      .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .

      x_n^(k+1) = b_n1x₁^(k+1) + b_n2x₂^(k+1) + b_n3x₃^(k+1) + … +  c_n .

Объединив приведение системы к виду, удобному для итераций и метод Зейделя в одну формулу, получим

      x_i^(k+1) = x_i^(k) – a_ii^–1(∑_j=1^i–1 a_ijx_j^(k+1) + ∑_j=1ⁿ a_ijx_i^(k) – b_i).

Тогда достаточным  условием сходимоти метода Зейделя  будет

      ∑_{j=1, j≠i} ⁿ | a_ij | < | a_ii |

(условие доминированния  диагонали).

      Метод Зейделя иногда называют также методом  Гаусса-Зейделя, процессом Либмана, методом последовательных замещений.

1.6. Сравнение прямых и итерационных методов

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

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

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

2.Практическая часть

2.1 Программа решения систем линейных уравнений по методу Гаусса

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

     Требуется решить систему линейных алгебраических уравнений с вещественными коэффициентами вида

      a₁₁x₁ + a₁₂x₂ + … + a_1nx_n = b₁ , 
 a₂₁x₂ + a₂₂x₂ + … + a_2nx_n = b₂ , 
 .   .   .   .   .   .   .   .   .   .   .   .   .

      a_n1x₁ + a_n2x₂ + … + a_nnx_n = b_n

для n ≤ 10 по методу Гаусса.

2.1.2. Тестовый пример

      3,2x₁ + 5,4x₂ + 4,2x₃ + 2,2x₄ = 2,6 ,

      2,1x₁ + 3,2x₂ + 3,1x₃ + 1,1x₄ = 4,8 ,

      1,2x₁ + 0,4x₂ – 0,8x₃ – 0,8x₄ = 3,6 ,

      4,7x₁ + 10,4x₂ + 9,7x₃ + 9,7x₄ = –8,4 ,

      x₁ = 5, x₂ = –4, x₃ = 3, x₄ = –2.

2.1.3. Описание алгоритма

     В данной программе реализован метод  Гаусса со схемой частичного выбора.

     В переменную n вводится  порядок матрицы системы. С помощью вспомогательной процедуры ReadSystem в двумерный массив a и одномерный массив b вводится c клавиатуры расширенная матрица системы, после чего оба массива и переменная n передаются функции Gauss. В функции Gauss для каждого k-го шага вычислений выполняется поиск максимального элемента в k-м столбце матрицы начиная с k-й строки. Номер строки, содержащей максимальный элемент сохраняется в переменной l. В том случае если максимальный элемент находится не в k-й строке, строки с номерами k и l меняются местами. Если же все эти элементы равны нулю, то происходит прекращение выполнения функции Gauss c результатом false. После выбора строки выполняется преобразование матрицы по методу Гаусса. Далее вычисляется решение системы и помещается в массив x. Полученное решение выводится на экран при помощи вспомогательной процедуры WriteX. 

      2.1.4. Листинг программы и результаты работы

Пример  №1:

Uses CRT;

Const

     maxn = 10;

Type

    Data = Real;

    Matrix = Array[1..maxn, 1..maxn] of Data;

    Vector = Array[1..maxn] of Data;

{ Процедура ввода  расширенной матрицы системы  }

Procedure ReadSystem(n: Integer; var a: Matrix; var b: Vector);

Var

   i, j, r: Integer;

Begin

     r := WhereY;

     GotoXY(2, r);

     Write('A');

     For i := 1 to n do begin

         GotoXY(i*6+2, r);

         Write(i);

         GotoXY(1, r+i+1);

         Write(i:2);

     end;

     GotoXY((n+1)*6+2, r);

     Write('b');

     For i := 1 to n do begin

         For j := 1 to n do begin

             GotoXY(j * 6 + 2, r + i + 1);

             Read(a[i, j]);

         end;

         GotoXY((n + 1) * 6 + 2, r + i + 1);

         Read(b[i]);

     end;

End;

{ Процедура вывода результатов }

Procedure WriteX(n :Integer; x: Vector);

Var

   i: Integer;

Begin

     For i := 1 to n do

         Writeln('x', i, ' = ', x[i]);

End; 
 

{ Функция, реализующая  метод Гаусса }

Function Gauss(n: Integer; a: Matrix; b: Vector; var x:Vector): Boolean;

Var

   i, j, k, l: Integer;

   q, m, t: Data;

Begin

   For k := 1 to n - 1 do begin

       { Ищем строку l с максимальным элементом в k-ом столбце}

         l := 0;

         m := 0;

         For i := k to n do

             If Abs(a[i, k]) > m then begin

                m := Abs(a[i, k]);

                l := i;

             end;

      { Если у всех строк от k до n элемент в k-м столбце нулевой,

                то система не имеет однозначного  решения }

         If l = 0 then begin

            Gauss := false;

            Exit;

         end;

   { Меняем местом l-ую строку с k-ой }

         If l <> k then begin

            For j := 1 to n do begin

                t := a[k, j];

                a[k, j] := a[l, j];

                a[l, j] := t;

            end;

            t := b[k];

            b[k] := b[l];

            b[l] := t;

         end;

  { Преобразуем матрицу }

         For i := k + 1 to n do begin

             q := a[i, k] / a[k, k];

             For j := 1 to n do

                 If j = k then

                    a[i, j] := 0

                 else

                      a[i, j] := a[i, j] - q * a[k, j];

                 b[i] := b[i] - q * b[k];

             end; 

     end;

     { Вычисляем решение }

     x[n] := b[n] / a[n, n];

     For i := n - 1 downto 1 do begin

         t := 0;

         For j := 1 to n-i do

             t := t + a[i, i + j] * x[i + j];

         x[i] := (1 / a[i, i]) * (b[i] - t);

     end;

     Gauss := true;

End; 

Var

    n, i: Integer;

    a: Matrix ;

    b, x: Vector;

Begin

      ClrScr;

      Writeln('Программа решения систем  линейных уравнений по методу  Гаусса');

      Writeln; 

      Writeln('Введите порядок матрицы системы (макс. 10)');

      Repeat

             Write('>');

             Read(n);

      Until (n > 0) and (n <= maxn);