Решение систем, линейных уравнений Опции

[Altair]

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

uses crt;

const
  nn = 10;
type
  Matrix = array[1..NN,1..NN+1] of real;

(* построчный ввод матрицы *)
Procedure ReadMatr(var A:Matrix; var n:word );
var
  i, j, m: word;
begin
  repeat
    write('Введите количество линейных уравн. в системе: '); readln(N)
  until (N>0) and (N<=NN);

  m:=n+1;
  For i:=1 to n do begin
    For j:=1 to m do begin
      write('A[',i,j,']= '); readln(A[i,j])
    end 
  end
end;

(* построчный вывод матрицы *)
Procedure PrintMatr(A:Matrix; n:word);
Var
  i, j, m: word;
begin
  m:=n+1;
  For i:=1 to n do begin
    For j:=1 to m do write(A[i,j],' ');
    writeln
  end
end;

procedure GaussM(a:matrix;n:word; var s:byte; var x:array of real);
var
  i, k, j: byte;
  m, t: real;
begin
  i:=1;
  s:=1;
  repeat
    j:=i+1;
    k:=i;
    m:=abs(a[i,i]);
    repeat
      if m<abs(a[j,i]) then begin
        m:=abs(a[j,i]);
        k:=j;
      end;
      j:=j+1
    until not(j<=n);

    if m<>0 then begin
      j:=i;
      repeat
        t:=a[i,j];
        a[i,j]:=a[k,j];
        a[k,j]:=t;
        j:=j+1
      until not(j<=n+1);
      k:=i+1;
      repeat
        t:=a[k,i]/a[i,i];
        a[k,i]:=0;
        j:=i+1;
        repeat
          a[k,j]:=a[k,j]-t*a[i,j];
          j:=j+1
        until not(j<=n+1);
        k:=k+1
      until not(k<=n);
    end
    else begin
      s:=0;
    end;
    i:=i+1
  until not((i<=n)and(s=1));

  if s=1 then begin
    i:=n;
    repeat
      x[i]:=a[i,n+1];
      j:=i+1;
      while j<=n do begin
        x[i]:=x[i]-a[i,j]*x[j];
        j:=j+1;
      end;
      x[i]:=x[i]/a[i,i];
      i:=i-1
    until not(i>=1);
  end;
end;

var
  b: array[0..nn] of real;
  a: Matrix;
  n, j: word;
  s: byte;
Begin
  readmatr(a,n);
  printmatr(a,n);
  writeln('press any key'); readkey;
  GaussM(a,n,s,b);
  for j:=1 to n do write (b[j],' ');
  writeln('press any key for exit ...'); readkey
end.

[volvo]

Решение систем линейных уравнений методом Крамера

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

{$N+}
{ Задаем порядок системы уравнений }
Const n = 3;

Type
  { Тип, описывающий матрицу системы (включая свободные члены !!!) }
  Equation = Array[1 .. n, 1 .. Succ(n)] Of Double;
  matrix =
    Array[1 .. n, 1 .. n] Of Double;

Const
  a: Equation =
   ((2, 1, 3, 9),
    (1, -2, 1, -2),
    (3, 2, 2, 7));

{ Процедура замены очередного столбца матрицы на свободные члены }
Procedure GetMatrix(wout: Integer; Var m: matrix);
  Var i, j: Integer;
  Begin
    For i := 1 To n Do
      For j := 1 To n Do
        If j <> wout Then m[i, j] := a[i, j]
        Else m[i, j] := a[i, Succ(n)]
  End;

{ Нахождение определителя }
function det(var p; const size: integer): double;

   function minusOne (n : integer): integer;
   begin
      minusOne := (1 - 2 * Byte (Odd (n)));
   end;
   function get_addr(i, j : integer;
            const n: integer) : integer;
   begin
      get_addr := pred (i) * n + j
   end;

   type vector = array[1 .. n * n] of double;

var
   my_p : vector absolute p;
   pp : ^vector;
   s : double;
   i, j, curr: integer;
begin
   s := 0.0;
   if size = 2 then
   begin
      det := my_p[1]*my_p[4] - my_p[2]*my_p[3]; exit
   end;

   for i := 1 to size do
   begin
      GetMem (pp, Sqr (pred (size)) * sizeof (double));
      curr := 1;
      for j := 1 to size do
         if j <> i then
         begin
            move (my_p[get_addr (j, 2, size)],
                  pp^[get_addr (curr, 1, pred (size))],
                  pred(size) * sizeof(double));
            inc (curr);
         end;

      s := s + minusOne (succ(i)) * my_p[get_addr (i, 1, size)] *
               det(pp^, pred (size));
      FreeMem (pp, Sqr (pred (size)) * sizeof (double))
   end;
   det := s
end;

Var
  i: Integer;
  mx: matrix;
  Determ: Double;
begin
  GetMatrix(Succ(n), mx);
  Determ := Det(mx, n);

  If Abs(Determ) < 1E-6 Then
    Writeln( 'Определитель исходной матрицы = 0' )
  Else
    For i := 1 To n Do
      Begin
        GetMatrix(i, mx);
        WriteLn( 'x(', i, ') = ', (Det(mx, n) / Determ):7:4 )
      End
end.