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

                      s:=s+r

                    end;

                  s:=a[i,i]-s;  if s<0 then l[i]:=-1 else l[i]:=1;

                  t[i,i]:=sqrt(abs(s));

                end;

            end;

          if i>j then t[i,j]:=0;

            if i<j then

              begin

                s:=0;i1:=i-1;

                for k:=1 to i1 do

                  begin

                    r:=t[k,i]*t[k,j];  if l[k]=-1 then r:=-r;  s:=s+r

                  end;

                t[i,j]:=(a[i,j]-s)/t[i,i];  if l[i]=-1 then t[i,j]:=-t[i,j];

              end

        end;{j}

    end;{i}

    for i:=1 to mmax do

      begin

        s:=0;i1:=i-1;

        for k:=1 to i1 do

          begin

            r:=t[k,i]*y[k];  if l[k]=-1 then r:=-r;  s:=s+r;

          end;

        y[i]:=(b[i]-s)/t[i,i];  if l[i]=-1 then y[i]:=-y[i];

      end;

     for i:=mmax downto 1 do

      begin

        s:=0;

        for k:=i+1 to mmax do

          s:=s+t[i,k]*x[k];

        x[i]:=(y[i]-s)/t[i,i];

      end;

end;

end.

Текст модуля (U_prog):

unit U_Prog;

Interface

uses UN_VVOD;

procedure Run_Prog(a:matrix; b:vector; n:integer; eps:real;

var x:vector; var stat:integer);

Implementation

var

   i, j  : integer;  hold_a, hold_b : vector;

procedure Run_Prog;

begin

  for i := 1 to n do

   if i = 1 then

     begin

       hold_a[i] := -A[1, 2] / A[1, 1];

       hold_b[i] := b[i] / A[1, 1];

     end

   else

     begin

       hold_a[i] := -A[i, i+1]/(A[i, i]+A[i, i-1]*hold_a[i-1]);

       hold_b[i] :=(-A[i, i-1]*hold_b[i-1] + b[i])/(A[i, i]+A[i,i-1]*hold_a[i-1]);

     end;

   x[n] := hold_b[n];

   for i := n-1 downto 1 do  x[i] := hold_a[i] * x[i+1] + hold_b[i];

   stat:=1;

end;  end.

Текст модуля (U_zeidel):

unit u_Zeidel;

interface

uses Un_VVod;

procedure Zeidel(A:matrix;B:vector;N:Integer;eps:real;

                        Var X:vector;var stat: integer);

implementation

Procedure Zeidel;

 Var C:matrix; D,Y:vector;  i,j,l:Integer; Delta:Real;

BEGIN

  For i:=1 to N do

    Begin

      d[i]:=0;

      For l:=1 to N do d[i]:=d[i]+a[l,i]*b[l];

      For j:=1 to N do

      Begin

        c[i,j]:=0;

        For l:=1 to N do c[i,j]:=c[i,j]+a[l,i]*a[l,j];

      End;

    End;

  For i:=1 to N do x[i]:=b[i];

  REPEAT

  Delta:=0;

  For i:=1 to n do

  Begin

    y[i]:=d[i]/c[i,i];

    For j:=1 to i-1 do y[i]:=y[i]-c[i,j]/c[i,i]*y[j];

    For j:=i+1 to N do y[i]:=y[i]-c[i,j]/c[i,i]*x[j];

    If Abs(y[i]-x[i])>Delta Then Delta:=Abs(y[i]-x[i]);

   x[i]:=y[i];

  End;

UNTIL Delta<=eps;

stat:=1;

END;  End.

Текст модуля (Un_vvod):

unit un_vvod;

interface

uses crt;

const n1=10;

type

     matrix=array[1..n1,1..n1] of real;

     vector=array[1..n1] of real;

var n:integer;

    a:matrix;

    b,x,e:vector;

    eps:real;

procedure VVOD(var n:integer;var eps:real;var a:matrix;var b:vector);

implementation

procedure VVOD;

var fn:string[15];

    f:text;

    i,j,k:integer;

    d,c:matrix;

begin

  writeln('Input data filename:');

  readln(fn);

  assign(f,fn);

  reset(f);

  readln(f,n,k,eps);

    {Reading Matrix D}

  for i:=1 to n do

    begin

      for j:=1 to n-1 do read(f,d[i,j]);

      readln(f,d[i,n]);

    end;

    {Reading Matrix C}

  for i:=1 to n do

    begin

      for j:=1 to n-1 do read(f,c[i,j]);

      readln(f,c[i,n]);

    end;

    {Reading Vector B}

  for i:=1 to n do read(f,b[i]);

     {Calculating Matrix A}

  for i:=1 to n do

    for j:=1 to n do a[i,j]:=d[i,j]+c[i,j]*k;

    close(f);

  end;

end.