Turbo Toolbox/Turbo_Toolbox.iso/sonderh1/linglei.pas

Pascal/Delphi Source File | 1986-10-12 | 6.3 KB | 247 lines

(* -------------------------------------------------------------------------
   Programm zum Loesen linearer Gleichungssysteme.
   Matrix und rechte Seite des Gleichungssystems werden aus einem
   Textfile eingelesen. Zeilen- und Spaltenzahl werden dabei
   automatisch bestimmt (Beispiel: WIDER.DAT).
   Zum Verfahren s.:
      Stoer J. , Heidelberger Taschenbuecher Nr.105:
      " Einfuehrung in die numerische Mathematik I ",
        S. 168 ff,
        Springer-Verlag , Berlin...

   Das Programm wurde auf einem ATARI 520ST+ mit
   Pascal ST plus  von CCD geschrieben.
   ------------------------------------------------------------------------- *)

PROGRAM lin_gleich (input, output);

CONST
  maxzeil = 60 ;
  maxspalt= 60 ;

TYPE
                         (* Manche Compiler koennen ueber dyn. Var. groessere
                              Felder erzeugen. Der von CCD allerdings nicht. *)
  ahi = array [1..maxzeil,1..maxspalt] OF REAL ;
  ptr = ^ahi ;

VAR
  zeil,                              (* Speichert die Anzahl der Gleichungen *)
  spalt : INTEGER;                   (*     "      "    "     "  Unbekannten *)

  a : ptr;              (* Feld fuer die Koeffizienten des Gleichungssystems *)
  d,                                  (* Hilfsfeld                           *)
  l : ARRAY [1..maxspalt] OF REAL;    (* Feld fuer das Ergebnis der Rechnung *)
  singmat : BOOLEAN;                  (* Matrix singulaer oder nicht ?!      *)

(* ------------------------------------------------------------------------- *)

PROCEDURE householder;

VAR
  sigma, s, sum, beta : REAL;
  i, j, k : INTEGER;

BEGIN
  Write ('Rechne');
  singmat := FALSE;
  j := 1;
  WHILE (j <= spalt) AND NOT(singmat) DO
  BEGIN
    Write ('.');                     (* Damit man sieht, dass etwas passiert *)
    sigma := 0;
    FOR i := j TO zeil DO
      sigma := sigma + a^[i,j] * a^[i,j];
    IF sigma <> 0 THEN
    BEGIN
      IF a^[j,j] < 0 THEN
      BEGIN
        s := Sqrt (sigma); d[j]:=s;
      END
      ELSE
      BEGIN
        s := -Sqrt (sigma); d[j]:=s;
      END;
      beta := 1/(s*a^[j,j]-sigma); a^[j,j] := a^[j,j] - s;
      FOR k := j+1 TO spalt+1 DO
      BEGIN
        sum := 0;
        FOR i := j TO zeil DO
          sum := sum + a^[i,j] * a^[i,k];
        sum := beta * sum ;
        FOR i := j TO zeil DO
          a^[i,k] := a^[i,k] + a^[i,j] * sum;
      END;
    END
    ELSE
    BEGIN
      WriteLn; WriteLn;
      WriteLn ('Singulaere Matrix');
      singmat := TRUE;
    END;
  j := Succ (j);
  END;
  WriteLn;
END; (* householder *)

(* ------------------------------------------------------------------------- *)

PROCEDURE ruecktransform;

VAR
   i, j : INTEGER;

BEGIN
  FOR i := spalt DOWNTO 1 DO
  BEGIN
    l[i] := a^[i,spalt+1];
    FOR j := spalt DOWNTO i+1 DO
    BEGIN
      l[i] := l[i] - a^[i,j] * l[j];
    END;
    l[i] := l[i]/d[i];
  END;
END; (* ruecktransform *)

(* ------------------------------------------------------------------------- *)

PROCEDURE ausgabe;

VAR
  i : INTEGER;

BEGIN
  WriteLn;
  IF zeil > spalt THEN
    Write ('Bestmoegliche ');
  Write ('Loesung des lin. Gleichungssystems: ');
  FOR i := 1 TO spalt DO
  BEGIN
    IF (i-1) MOD 3 = 0 THEN
      WriteLn;
    Write ('a[',pred(i):2,']',' = ', l[i]:14,'    ');
  END;
  WriteLn;
  WriteLn;
END; (* ausgabe *)

(* -------------------------------------------------------------------------
   Diese Prozedur wird im Moment nicht benutzt! Mit ihr kann die Matrix
   waehrend der Rechnung auf dem Bildschirm ausgegeben werden.               *)

PROCEDURE schrmatrix;

VAR
  i, j : INTEGER;

BEGIN
  WriteLn;
  WriteLn ('HOUSEHOLDERMATRIX :');
  FOR i := 1 TO spalt DO
  BEGIN
    FOR j := 1 TO spalt+1 DO
      Write (a^[i,j]:6,' ');
    WriteLn;
  END;
END; (* schrmatrix *)

(* ------------------------------------------------------------------------- *)

PROCEDURE pagen;           (* Bildschirm loeschen und Cursor nach links oben *)

BEGIN
  Write (Chr(27));
  Write ('E');
END; (* pagen *)

(* -------------------------------------------------------------------------
   Diese Prozedur kann auch so geschrieben werden, dass kein zweimaliges
   Lesen des Datenfiles noetig ist. Dies wird dann aber eher unuebersicht-
   lich.
   Impl.-Hinweis: Reset (filvar, 'filename')
                  oeffnet die Datei 'filename' zum lesen.                    *)

PROCEDURE einlesen;

VAR
  name    : STRING [25];
  tex     : TEXT;
  n, i, j,
  numcount: INTEGER;
  buff    : STRING [240];
  c1, c2  : CHAR;

BEGIN
  pagen;
  Write ('Name des Datenfiles: ');
  ReadLn (name);
  WriteLn;
  Write ('Lese Daten...');
  Reset (tex, name);
  REPEAT                                                     (* BEGIN suchen *)
    ReadLn (tex, buff);
  UNTIL Pos ('BEGIN', buff) = 1;
  spalt := 0;                  (* Spaltenzahl + 1 aus der 1. Zeile bestimmen *)
  c1 := ' ';
  ReadLn (tex, buff);
  FOR i := 1 TO Length (buff) DO
  BEGIN
    c2 := buff[i];
    IF ((c1=' ') AND (c2<>' ')) THEN
      spalt := Succ (spalt);
    c1 := c2;
  END;
  numcount := 0;
  REPEAT
    ReadLn (tex, buff);
    numcount := Succ (numcount);                    (* dabei  Zeilen zaehlen *)
  UNTIL Pos ('END', buff) = 1;
  Reset (tex, name);                                (* Zurueck an den Anfang *)
  REPEAT
    ReadLn (tex, buff);
  UNTIL Pos ('BEGIN', buff) = 1;
  i := 0; n := 0;
  REPEAT                                           (* Koeffizienten einlesen *)
    n := Succ(n);
    FOR i := 1 TO spalt DO
      Read (tex, a^[n,i]);
    Write ('.');
  UNTIL n = numcount;
  zeil := numcount;                                  (* Zeilenzahl festlegen *)
  spalt := Pred (spalt);                          (* Spaltenzahl korrigieren *)
  WriteLn; WriteLn;
END; (* einlesen *)

(* ------------------------------------------------------------------------- *)

PROCEDURE ende;

VAR
  anything : CHAR;

BEGIN
  Write ('Mit <RETURN> zurueck zum Desktop');
  Read (anything);
END; (* ende *)

(* ------------------------------------------------------------------------- *)

BEGIN (* lin_gleich *)
  New (a);
  einlesen;
  householder;
  IF NOT(singmat) THEN
  BEGIN
    ruecktransform;
    ausgabe;
  END;
  ende;
end.