Наиболее известным и популярным точным способом решения систем линейных алгебраических уравнений (СЛАУ) является метод Гаусса. Этот метод заключается в последовательном исключении неизвестных. Пусть в системе уравнений

первый элемент a₁₁⁽⁰⁾ не равен 0. Назовем его ведущим элементом первой строки. Поделим все элементы этой строки на a₁₁⁽⁰⁾ и исключим x₁ из всех последующих строк, начиная со второй, путем вычитания первой (преобразованной), умноженной на коэффициент при x₁ в соответствующей строке. Получим

Если a₂₂⁽¹⁾, то, продолжая аналогичное исключение, приходим к системе уравнений с верхней треугольной матрицей

Из нее в обратном порядке находим все значения x_i:

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

{ Написать программу решения СЛАУ Ax=b методом Гаусса. }

program SLAU;
const size=30;
type matrix=array [1..size,1..size] of real;
type vector=array [1..size] of real;

procedure GetDimension (var n:integer; s:string);
begin
 repeat
  writeln (s);
  readln (n);
  if (n<2) or (n>size) then writeln ('Должно быть число от 2 до ', size);
 until (n>1) and (n<size);
end;

procedure GetMatrix (n,m:integer; var a:matrix);
var i,j:integer;
begin
 for i:=1 to n do begin
  writeln ('Строка ',i);
  for j:=1 to m do read (a[i,j]);
 end;
end;

procedure PutMatrix (n,m:integer; var a:matrix);
var i,j:integer;
begin
 for i:=1 to n do begin
  writeln;
  for j:=1 to m do write (a[i,j]:10:3);
 end;
end;

procedure GetVector (n:integer; var a:vector);
var i:integer;
begin
 for i:=1 to n do begin
  write (i,' элемент=');
  readln (a[i]);
 end;
end;

procedure PutVector (n:integer; var a:vector);
var i:integer;
begin
 writeln;
 for i:=1 to n do writeln (a[i]:10:3);
end;

procedure DirectGauss (n:integer; var a:matrix; var b:vector);
var i,j,k:integer;
    buf:real;
begin
 for i:=1 to n-1 do
 for j:=i+1 to n do begin
  buf:=a[i,i]/a[j,i];
  for k:=1 to n do a[j,k]:=a[j,k]*buf-a[i,k];
  b[j]:=b[j]*buf-b[i];
 end;
end;

procedure ObratGauss (n:integer; var a:matrix; var x,b:vector);
var i,j,k:integer;
    buf:real;
begin
 x[n]:=b[n]/a[n,n];
 for i:=n-1 downto 1 do begin
  buf:=0;
  for j:=i+1 to n do buf:=buf+a[i,j]*x[j];
  x[i]:=(b[i]-buf)/a[i,i];
 end;
end;

procedure Mvmul (n,m:integer; var a:matrix; var x,b:vector);
var i,j:integer;
begin
 for i:=1 to n do begin
  b[i]:=0;
  for j:=1 to m do b[i]:=b[i]+a[i,j]*x[j];
 end;
end;

var a,a2:matrix;
    x,b:vector;
    n,m:integer;

begin
 GetDimension (n,'Введите размерность квадратной матрицы');
 m:=n;
 GetMatrix (n,m,a);
 a2:=a;
 writeln ('Ввод правой части:');
 GetVector (n,b);
 DirectGauss (n,a,b);
 ObratGauss (n,a,x,b);
 write ('Решение:');
 PutVector (m,x);
 write ('Проверка A*x:');
 Mvmul (n,m,A2,x,b);
 PutVector (n,b);
 Reset (Input); Readln;
end.

Аналогичная программа на C++ выглядит следующим образом:

/*
 Написать программу для решения системы линейных алгебраческих уравнений (СЛАУ)
 методом Гаусса
*/
#include <iostream.h>
#include <conio.h>
#include <stdlib.h>
 
void main () {
 const int n=3; //Размерность системы
 float a[n][n+1],a0[n][n+1];
  //A - расширенная матрица системы, A0 - ее копия для проверки решения
 int i,j,k;
 float buf,x[n],b[n];
 clrscr();
 randomize ();
 cout << "Матрица и вектор правой части";
 for (i=0; i<n; i++) {
  cout << endl;
  for (j=0; j<n+1; j++) {
   a0[i][j]=a[i][j]=1+random(100)/25.; 
   //Нулевых элементов нет, деления на ноль не будет
   cout << a[i][j] << " ";
  }
 }
 
 //Прямой ход метода Гаусса
 for (i=0;i<n-1;i++)
 for (j=i+1;j<n;j++) {
  buf=a[i][i]/a[j][i];
  for (k=0;k<=n;k++) a[j][k]=a[j][k]*buf-a[i][k];
 }
 //Обратный ход метода Гаусса
 x[n-1]=a[n-1][n]/a[n-1][n-1];
 for (i=n-2;i>=0;i--) {
  buf=0;
  for (j=i+1;j<n;j++) buf+=a[i][j]*x[j];
  x[i]=(a[i][n]-buf)/a[i][i];
 }
 
 cout << endl << "Решение" << endl;
 for (i=0; i<n; i++) cout << x[i] << " ";
 
 cout << endl << "Проверка" << endl;
 for (i=0; i<n; i++) {
  b[i]=0;
  for (j=0; j<n; j++) b[i]+=a0[i][j]*x[j];
  cout << b[i] << " ";
 }
 
}

Здесь матрица и вектор правой части генерируются случайным образом из чисел в диапазоне от 1 до 5:

1+random(100)/25.

Естественно, диапазон можно поменять.

Лучшее решение на Паскале есть вот здесь.

гостевая; E-mail