[angvelem]

Примерно так:

Код:

{*** Решение уравнений методом Гаусса - Жордана ***}
unit GausJordan;

interface

type
  TInMatr		= array[0..19] of array[0..19] of Double;
  TOutMatr		= array[0..19] of Double;

const
  GAUSS_ACCURACY	= 0.0000001;
  GAUSS_OK		= 0;
  GAUSS_NOSOL		= 1;
  GAUSS_MANYSOL		= 2;

function gauss(a : TInMatr; n, u : Integer; var x : TOutMatr) : Integer;

implementation

{********************************************* 
 * Функция решения систем линейных уравнений * 
 * методом Гаусса - Жордана                  * 
 * (C) 2002 Восков Алексей                   * 
 * версия 2.1                                * 
 *********************************************}
// ВХОДНЫЕ ДАННЫЕ
{  a[20][20] - Матрица для хранения системы a[y][x] 
           последний столбец - для хранения св. члена 
           св. член - в правой части 
   n - число неизвестных    
   u - число уравнений    
   ВЫХОДНЫЕ ДАННЫЕ    
   x[20] - массив для хранения корней системы      
   в случае нормального выполнения задачи функция возвращает 0, 
   в случае неразрешимости 1   в случае бесконечного числа решений 2 }

function gauss(a : TInMatr; n, u : Integer; var x : TOutMatr) : Integer;
var 
  i, j, k : Integer;		// Счетчики циклов
  sn      : Integer;		// Номер строки
  d       : Double;		// Коэффициент домножения или модуль наиб. эл.
begin  
  Result := GAUSS_OK;		// Нормальное завершение работы
  
  //*** Проверка u и n ***
  if n > u then
  begin
    Result := GAUSS_MANYSOL;
    Exit;
  end;
  
  //*** Приведение к диагональному виду ***
  for j := 0 to n - 1 do
  begin         
    // а) поиск строки с наибольшим по модулю элементом
    d  := Abs(a[j][j]);
    sn := j;

    for i := j to u - 1 do
      if Abs(a[i][j]) > d then
      begin
        d  := Abs(a[i][j]);
        sn := i;
      end;
       
    // б) перенос строки на надлежащее место
    for k := 0 to n do
    begin
      d := a[sn][k];
      a[sn][k] := a[j][k];
      a[j][k]  := d;
    end;
     
    // в) деление ведущего ур-я на главный элемент
    d := a[j][j];
     
    if d > 0 then
      for k := 0 to n do
        a[j][k] := a[j][k] / d
    else
      for k := 0 to n do
        a[j][k] := 0;
       
    // г) вычитание данной строки из всех остальных
    //    с домножением на коэффициент
    for i := 0 to u - 1 do
    begin
      if i = j then
        Continue;			// Не трогаем вычит. строку
      d := -a[i][j];

      for k := 0 to n do		// Вычитание
        a[i][k] := a[i][k] + a[j][k] * d;
    end;
  end;

  //*** Вычисление корней ***
  // а) проверка системы на разрешимость
  if u > n then
  begin
    for i := n to u - 1 do
    begin
      k := 0;
      for j := 0 to n - 1 do
        if Abs(a[i][j]) > GAUSS_ACCURACY then
          k := 1;
 
      if (k = 0) and (Abs(a[i][n]) > GAUSS_ACCURACY) then
      begin
        Result := GAUSS_NOSOL;
        Exit;
      end;
    end;
  end;
   
  // б) поиск корней
  for i := 0 to n - 1 do
  begin
    x[i] := -a[i][n];

    if a[i][i] <> 1 then		// Обработка ошибок
    begin
      if x[i] > 0 then
      begin
        Result := GAUSS_NOSOL;		// Решений нет
        Exit;
      end
      else
      begin
        Result := GAUSS_MANYSOL;	// Бесконечно много решений
        Exit;
      end;
    end;

    if Abs(x[i]) < GAUSS_ACCURACY then
      x[i] := 0;			// Обнуление слишком малых знач.
  end;
end;

end.