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

[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.

[hiv]

Решение системы линейных уравнений методом Холецкого

Этот метод менее подвержен накопительным ошибкам в сравнении с методом Гаусса. Поэтому им можно решать системы из многих десятков уравнений!

Код

program HOLETSKY;

{$M 65520, 65520, 65520}

const max_eq = 50;
type
   RealType=double;

   MCvector = array [0..max_eq] of RealType; {вектор}

   MCmatr = array [0..max_eq,0..max_eq] of RealType; {запись типа матрица}

function maxval(x,y:realtype):realtype;
begin
    if x<y then maxval:=y
    else maxval:=x;
end;

procedure Holetskiy(N:integer;A0:MCmatr; B0:MCvector; var X:MCvector; var error:RealType);
var
    i,j,m,k,p :byte;
    s,e,max   :RealType;
    A :MCmatr;
begin
    if (N>max_eq)or(N<2) then
    begin
         writeln('Ошибка в передаваемых данных.');
         halt(1);
    end;

    A:=A0;
    for i:=0 to N-1 do A[i,N]:=B0[i];

    for i:=0 to N-1 do       {выбор главного элемента по диагонали}
    begin
         max:=0;
         for j:=i to N-1 do
             if max<abs(A[i,j]) then
             begin
                  max:=abs(A[i,j]);
                  p:=j;
             end;
         if p<>i then
         begin
              for j:=0 to N do A[N,j]:=A[i,j];
              for j:=0 to N do A[i,j]:=A[p,j];
              for j:=0 to N do A[p,j]:=A[N,j];
         end;
    end;

    {метoд  хoлецкoгo}
    for j:=1 to N do A[0,j]:=A[0,j]/A[0,0];
    for m:=1 to N-1 do
    begin
         for i:=m to N-1 do
         begin
              s:=0;
              for k:=0 to m-1 do s:=s+A[i,k]*A[k,m];
              A[i,m]:=A[i,m]-s
         end;
         for j:=m+1 to N do
         begin
              s:=0;
              for k:=0 to m-1 do s:=s+A[m,k]*A[k,j];
              A[m,j]:=(A[m,j]-s)/A[m,m]
         end;
    end;
    X[N-1]:=A[N-1,N];
    for i:=N-2 downto 0 do
    begin
         s:=0;
         for k:=i+1 to N-1 do s:=s+A[i,k]*X[k];
         X[i]:=A[i,N]-s
    end;

    {расчет ошибки}
    error:=0;
    for i:=0 to N-1 do
    begin
         s:=0;
         for j:=0 to N-1 do
             if i=j then s:=s+(A0[i,j]+1)*X[j]
             else s:=s+A0[i,j]*X[j];
         s:=s-B0[i];
         error:=maxval(error,abs(X[i]-s));
    end;
end;

var A   :MCmatr;
   B,X :MCvector;
   err :RealType;
   n,i,j   :integer;

begin {of program}
    writeln('Решение системы линейных уравнений метoдoм Хoлецкoгo');
    writeln('с выбoрoм главнoгo элемента.');
    write('Введите числo уравнений = ');read(N);
    for i:=0 to N-1 do
    begin
         writeln('Введите кoэффициенты ',i:0,'-го уравнения:');
         for j:=0 to N-1 do
         begin
           write('A[',i:0,',',j:0,'] = '); readln(A[i,j]);
         end;
         write('B[',i:0,'] = '); readln(B[i]);
    end;

    Holetskiy(N,A,B,X,err);

    writeln('Решение уравнения:');
    for i:=0 to N-1 do  writeln('X[',i:0,'] = ',X[i]);
    writeln('Точность резутьтата =',err);
end.