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

uses crt;
const raz=4;
type matrix=array[1..raz,1..raz+1] of real;
var matr:matrix;
    x:array[1..raz] of real;
    m,p:integer;
    s,b:real;
    matr_a,matr_b,matr_c,matr_d:matrix;
 
procedure Vivod(matr:matrix);
begin
  for m:=1 to raz do begin
      for p:=1 to (raz+1) do begin
          write(matr[m,p]:8:1);
      end;
      writeln;
  end;
end;
 
procedure Gauss;
var i,j,k,l:byte;
    s,b,a_matr:real;
    matr_a:matrix;
begin
  matr_a:=matr;
  for j:=1 to raz do begin
      a_matr:=matr_a[j,j];
      for i:=1 to (raz+1) do
      matr_a[j,i]:=matr_a[j,i]/a_matr;
      if j<>raz then
        for k:=1 to (raz-j) do begin
          b:=matr_a[j+k,j];
          for i:=1 to (raz+1) do
            matr_a[j+k,i]:=matr_a[j+k,i]-matr_a[j,i]*b;
        end;
    for l:=1 to raz do
      if matr_a[l,l]=0 then begin
        writeln('Na glavnoy diagonali est nulevie elementi !');halt;
      end;
  end;
 
  for k:=raz downto 1 do begin
    s:=0;
    for i:=1 to raz do begin
      s:=s-matr_a[k,i]*x[i];
    end;
    x[k]:=s+matr_a[k,raz+1];
  end;
  writeln;
  writeln('Metod Gaussa:');
  for i:=1 to raz do
    if (i mod 2)<>0 then
        write('X',i,'=',x[i]:7:3,'   ')
    else
        writeln('X',i,'=',x[i]:7:3,'   ');
end;
 
procedure Zeydel;
  var i,j,k,l:integer;
      s,e:real;
begin
 
  for i:=1 to raz do x[i]:=0;
  matr_c:=matr;
  e:=1e-2;
  m:=100;
 
  for k:=1 to m do begin
    l:=k;
    for i:=1 to raz do begin
      s:=matr_c[i,raz+1];
 
      for j:=1 to raz do
        s:=s-matr_c[i,j]*x[j];
      s:=s/matr_c[i,i];
      x[i]:=x[i]+s;
      if abs(s)>e then l:=0
    end;
    if l<>0 then Break;
  end;
  writeln;
  writeln('Metod Zeydela');
  for i:=1 to raz do
    if (i mod 2)<>0 then
        write('X',i,'=',x[i]:7:3,'   ')
    else
        writeln('X',i,'=',x[i]:7:3,'   ');
end;
 
procedure Progon;
var i,j:integer;
    alfa,betta,s1,s2,s3:array[1..raz] of real;
    xi1,xin:real;
begin
  matr_d:=matr;
 
  for i:=1 to raz do x[i]:=0;
 
  alfa[1]:=-matr_d[1,2]/matr_d[1,1];
  betta[1]:=matr_d[1,raz+1]/matr_d[1,1];
  xi1:=-matr_d[1,2]/matr_d[1,1];
  xin:=-matr_d[raz,raz-1]/matr_d[raz,raz];
 
  for i:=2 to raz do begin
    for j:=1 to raz-1 do begin
      s1[j]:=matr_d[j+1,j];
      s2[j]:=matr_d[j,j+1];
      s3[j]:=matr_d[j,j];
    end;
    alfa[i]:=-matr_d[i,i+1]/(matr_d[i,i]+alfa[i-1]*matr_d[i+1,i]);
    betta[i]:=(-betta[i-1]*matr_d[i+1,i]+matr_d[i,raz+1])/
              (matr_d[i,i]+alfa[i-1]*matr_d[i+1,i]);
  end;
  x[raz]:=(-matr_d[raz,raz-1]*betta[raz-1]+matr_d[raz,raz+1])/
             (matr_d[raz,raz]+alfa[raz-1]*matr_d[raz,raz-1]);
  for j:=raz-1 downto 1 do begin
    x[j]:=alfa[j]*x[j+1]+betta[j];
  end;
  if (abs(xi1)<=1) and (abs(xin)<1) and
     (s1[j]<>0) and (s2[j]<>0) and
     (abs(s3[j])>=(abs(s1[j])+abs(s2[j]))) then
  begin
    writeln;
    writeln('Metod Progonki:');
  for i:=1 to raz do
    if (i mod 2)<>0 then
        write('X',i,'=',x[i]:7:3,'   ')
    else
        writeln('X',i,'=',x[i]:7:3,'   ');
 
  end
  else writeln('Progonka ne ustoychiva');
end;
 
procedure LU_razl;
var i,j,k:byte;
begin
  matr_b:=matr;
{ L-U }
  for j:=1 to raz do
    for i:=1 to raz do
      if i<=j then begin
        s:=0;
        if i>1 then
          for k:=1 to (i-1) do
            s:=s+matr_b[j,k]*matr_b[k,i];
        matr_b[j,i]:=matr_b[j,i]-s;
      end
      else begin
        s:=0;
        if j>1 then
          for k:=1 to (j-1) do
            s:=s+matr_b[j,k]*matr_b[k,i];
        matr_b[j,i]:=(matr_b[j,i]-s)/matr_b[j,j];
      end;
{ Z }
  for j:=1 to raz do begin
    s:=0;
    if j>1 then
      for i:=1 to (j-1) do
        s:=s+matr_b[j,i]*matr_b[i,raz+1];
    matr_b[j,raz+1]:=(matr_b[j,raz+1]-s)/matr_b[j,j];
  end;
{ X }
  for j:=raz downto 1 do begin
    s:=0;
    if j<raz then
      for i:=(j+1) to raz do
        s:=s+matr_b[j,i]*matr_b[i,raz+1];
    matr_b[j,raz+1]:=matr_b[j,raz+1]-s;
  end;
  writeln;
  writeln('Metod LU razlojeniya:');
  for i:=1 to raz do
    if (i mod 2)<>0 then
        write('X',i,'=',matr_b[i,raz+1]:7:3,'   ')
    else
        writeln('X',i,'=',matr_b[i,raz+1]:7:3,'   ');
end;
 
Begin
  clrscr;
  matr[1,1]:=-3;matr[1,2]:=1;matr[1,3]:=0;matr[1,4]:=0;matr[1,5]:=3.9;
  matr[2,1]:=-1;matr[2,2]:=3.6;matr[2,3]:=4;matr[2,4]:=0;matr[2,5]:=4.2;
  matr[3,1]:=0;matr[3,2]:=-1.6;matr[3,3]:=1.5;matr[3,4]:=5;matr[3,5]:=1;
  matr[4,1]:=0;matr[4,2]:=0;matr[4,3]:=-2;matr[4,4]:=-15;matr[4,5]:=2;
 
  Vivod(matr);
    writeln('Press any key');
    readkey;
  Gauss;
  LU_razl;
  Zeydel;
  Progon;
 
  readkey;
end.