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

Введение

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

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

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

Известны примеры решенных в  последние годы задач, где число неизвестных достигало сотен тысяч. Естественно, это было бы невозможно, если бы соответствующие матрицы не являлись разреженными (матрица системы из 100 тыс. уравнений в формате двойной точности заняла бы около 75 Гбайт).

1. Теоретическая часть

1.1. Метод Гаусса

Одним из самых распространенных методов решения систем линейных уравнений является метод Гаусса. Этот метод (который также называют методом последовательного исключения неизвестных) известен в различных вариантах уже более 2000 лет.

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

1. Метод Гаусса с выбором главного элемента по столбцу

(схема частичного  выбора). Описание метода. На k-м шаге прямого хода коэффициенты уравнений системы с номерами i = k + 1, …, n преобразуются по формулам

a_ij^(k) = a_ij^(k–1) − q_ika_kj , b_i^(k) = b_i^(k–1) − q_ikb_k^(k–1) , i = k + 1, …, n.

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

В методе Гаусса с выбором  главного элементоа по столбцу гарантируется, что |q_ik| ≤ 1 для всех k = 1, 2, …, n – 1 и i = k + 1, …, n. Отличие этого варианта метода Гаусса от схемы единственного деления заключается в том, что на k-м шаге исключения в качестве главного элемента выбирают максимальный по модулю коэффициент a_ikk при неизвестной x_k в уравнениях с номерами i = k + 1, …, n. Затем соответствующее выбранному коэффициенту уравнение с номером i_k меняют местами с k-м уравнением системы для того, чтобы главный элемент занял место коэффициента a_kk^(k-1). После этой перестановки исключение неизвестного x_k производят, как в схеме единственного деления.

 

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 по методу Гаусса c выбором ведущего элемента по столбцу.

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. Описание  алгоритма

 

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

{******************************************************************************

 

Решение систем линейных алгебраических уравнений (СЛАУ)

Метод Гаусса. Выбор ведущего элемента по столбцу v. 1.0

*************************************************************

 

uses crt;

 

const

   MAX_SIZE=20;                          {максимальный размер матрицы}

 

type

   TMatrix = array [1..MAX_SIZE] of array [1..MAX_SIZE+1] of double;

   TVector = array [1..MAX_SIZE] of double;

var

   a: TMatrix;

   x: TVector;

   n, i, j: integer;

   data_file: text;

   result_file: text;

 

{обмен строк матрицы}

procedure swapRow(var a: TMatrix;

                  n: integer;

                  i1: integer;

                  i2: integer);

var

   j: integer;

   t: double;

begin

   for j:=1 to n+1 do begin

      t := a[i1,j];

      a[i1,j] := a[i2,j];

      a[i2,j] := t;

   end

end;

 

 

{решение СЛАУ методом Гаусса}

procedure linearSolveGaussColLead(a: TMatrix;     {расширенная матрица коэфф.}

                                  var x: TVector; {результат}

                                  n: integer);

var

   i, j, k: integer;

   max, lead, a_div_lead, sum, eps: double;

   max_i : integer;

 

begin

   eps := 0.000001;

 

   {прямой ход}

   for k:=1 to n do begin

      {поиск максимального по модулю элемента в столбце}

      max := 0;

      max_i := -1;

 

      for i:=k to n do begin

if(abs(a[i,k]) > max) then begin

    max := abs(a[i,k]);

    max_i := i;

         end;

      end;

 

      if((max_i = -1) or (abs(a[max_i,k])<eps)) then begin

 writeln('Ошибка: определитель равен нулю');

      end

      else begin

if(max_i <> k) then begin

    {обмен строк}

    swapRow(a,n,max_i,k);

         end;

      end;

 

      {обнуление элементов}

      lead := a[k,k];

 

      for j:=k to n+1 do begin

a[k,j] := a[k,j]/lead;

      end;

 

      for i:=k+1 to n do begin

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

 

for j:=k to n+1 do begin

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

end;

      end;

   end;

 

   {обратный ход}

   x[n] := a[n][n+1];

 

   for k:=n-1 downto 1 do begin

      sum := a[k][n+1];

 

      for j:=k+1 to n do begin

sum := sum-a[k,j]*x[j];

      end;

 

      x[k] := sum;

   end;

end;

 

 

begin

   clrscr;

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

   writeln('Метод Гаусса. Выбор ведущего элемента по столбцу');

   writeln;

 

   assign(data_file, 'a.txt');

   reset(data_file);

   read(data_file, n);

 

   for i:=1 to n do begin

      for j:=1 to n+1 do begin

         read(data_file, a[i,j]);

      end;

   end;

 

   close(data_file);

 

   writeln('Расширенная матрица системы:');

   for i:=1 to n do begin

      for j:=1 to n+1 do begin

         write(a[i,j]:14:6);

      end;

      writeln;

   end;

   writeln;

 

   linearSolveGaussColLead(a,x,n);

 

   writeln('Результат:');

   for i:=1 to n do begin

      writeln(x[i]:14:6);

   end;

   writeln;

 

   assign(result_file, 'result.txt');

   rewrite(result_file);

 

   for i:=1 to n do begin

      writeln(result_file, x[i]:14:6);

   end;

 

   close(result_file);

 

   readkey;

end.