home *** CD-ROM | disk | FTP | other *** search
/ Turbo Toolbox / Turbo_Toolbox.iso / sonderh1 / linglei.pas < prev    next >
Pascal/Delphi Source File  |  1986-10-12  |  6KB  |  247 lines
  1. (* -------------------------------------------------------------------------
  2.    Programm zum Loesen linearer Gleichungssysteme.
  3.    Matrix und rechte Seite des Gleichungssystems werden aus einem
  4.    Textfile eingelesen. Zeilen- und Spaltenzahl werden dabei
  5.    automatisch bestimmt (Beispiel: WIDER.DAT).
  6.    Zum Verfahren s.:
  7.       Stoer J. , Heidelberger Taschenbuecher Nr.105:
  8.       " Einfuehrung in die numerische Mathematik I ",
  9.         S. 168 ff,
  10.         Springer-Verlag , Berlin...
  11.  
  12.    Das Programm wurde auf einem ATARI 520ST+ mit
  13.    Pascal ST plus  von CCD geschrieben.
  14.    ------------------------------------------------------------------------- *)
  15.  
  16. PROGRAM lin_gleich (input, output);
  17.  
  18. CONST
  19.   maxzeil = 60 ;
  20.   maxspalt= 60 ;
  21.  
  22. TYPE
  23.                          (* Manche Compiler koennen ueber dyn. Var. groessere
  24.                               Felder erzeugen. Der von CCD allerdings nicht. *)
  25.   ahi = array [1..maxzeil,1..maxspalt] OF REAL ;
  26.   ptr = ^ahi ;
  27.  
  28. VAR
  29.   zeil,                              (* Speichert die Anzahl der Gleichungen *)
  30.   spalt : INTEGER;                   (*     "      "    "     "  Unbekannten *)
  31.  
  32.   a : ptr;              (* Feld fuer die Koeffizienten des Gleichungssystems *)
  33.   d,                                  (* Hilfsfeld                           *)
  34.   l : ARRAY [1..maxspalt] OF REAL;    (* Feld fuer das Ergebnis der Rechnung *)
  35.   singmat : BOOLEAN;                  (* Matrix singulaer oder nicht ?!      *)
  36.  
  37. (* ------------------------------------------------------------------------- *)
  38.  
  39. PROCEDURE householder;
  40.  
  41. VAR
  42.   sigma, s, sum, beta : REAL;
  43.   i, j, k : INTEGER;
  44.  
  45. BEGIN
  46.   Write ('Rechne');
  47.   singmat := FALSE;
  48.   j := 1;
  49.   WHILE (j <= spalt) AND NOT(singmat) DO
  50.   BEGIN
  51.     Write ('.');                     (* Damit man sieht, dass etwas passiert *)
  52.     sigma := 0;
  53.     FOR i := j TO zeil DO
  54.       sigma := sigma + a^[i,j] * a^[i,j];
  55.     IF sigma <> 0 THEN
  56.     BEGIN
  57.       IF a^[j,j] < 0 THEN
  58.       BEGIN
  59.         s := Sqrt (sigma); d[j]:=s;
  60.       END
  61.       ELSE
  62.       BEGIN
  63.         s := -Sqrt (sigma); d[j]:=s;
  64.       END;
  65.       beta := 1/(s*a^[j,j]-sigma); a^[j,j] := a^[j,j] - s;
  66.       FOR k := j+1 TO spalt+1 DO
  67.       BEGIN
  68.         sum := 0;
  69.         FOR i := j TO zeil DO
  70.           sum := sum + a^[i,j] * a^[i,k];
  71.         sum := beta * sum ;
  72.         FOR i := j TO zeil DO
  73.           a^[i,k] := a^[i,k] + a^[i,j] * sum;
  74.       END;
  75.     END
  76.     ELSE
  77.     BEGIN
  78.       WriteLn; WriteLn;
  79.       WriteLn ('Singulaere Matrix');
  80.       singmat := TRUE;
  81.     END;
  82.   j := Succ (j);
  83.   END;
  84.   WriteLn;
  85. END; (* householder *)
  86.  
  87. (* ------------------------------------------------------------------------- *)
  88.  
  89. PROCEDURE ruecktransform;
  90.  
  91. VAR
  92.    i, j : INTEGER;
  93.  
  94. BEGIN
  95.   FOR i := spalt DOWNTO 1 DO
  96.   BEGIN
  97.     l[i] := a^[i,spalt+1];
  98.     FOR j := spalt DOWNTO i+1 DO
  99.     BEGIN
  100.       l[i] := l[i] - a^[i,j] * l[j];
  101.     END;
  102.     l[i] := l[i]/d[i];
  103.   END;
  104. END; (* ruecktransform *)
  105.  
  106. (* ------------------------------------------------------------------------- *)
  107.  
  108. PROCEDURE ausgabe;
  109.  
  110. VAR
  111.   i : INTEGER;
  112.  
  113. BEGIN
  114.   WriteLn;
  115.   IF zeil > spalt THEN
  116.     Write ('Bestmoegliche ');
  117.   Write ('Loesung des lin. Gleichungssystems: ');
  118.   FOR i := 1 TO spalt DO
  119.   BEGIN
  120.     IF (i-1) MOD 3 = 0 THEN
  121.       WriteLn;
  122.     Write ('a[',pred(i):2,']',' = ', l[i]:14,'    ');
  123.   END;
  124.   WriteLn;
  125.   WriteLn;
  126. END; (* ausgabe *)
  127.  
  128. (* -------------------------------------------------------------------------
  129.    Diese Prozedur wird im Moment nicht benutzt! Mit ihr kann die Matrix
  130.    waehrend der Rechnung auf dem Bildschirm ausgegeben werden.               *)
  131.  
  132. PROCEDURE schrmatrix;
  133.  
  134. VAR
  135.   i, j : INTEGER;
  136.  
  137. BEGIN
  138.   WriteLn;
  139.   WriteLn ('HOUSEHOLDERMATRIX :');
  140.   FOR i := 1 TO spalt DO
  141.   BEGIN
  142.     FOR j := 1 TO spalt+1 DO
  143.       Write (a^[i,j]:6,' ');
  144.     WriteLn;
  145.   END;
  146. END; (* schrmatrix *)
  147.  
  148. (* ------------------------------------------------------------------------- *)
  149.  
  150. PROCEDURE pagen;           (* Bildschirm loeschen und Cursor nach links oben *)
  151.  
  152. BEGIN
  153.   Write (Chr(27));
  154.   Write ('E');
  155. END; (* pagen *)
  156.  
  157. (* -------------------------------------------------------------------------
  158.    Diese Prozedur kann auch so geschrieben werden, dass kein zweimaliges
  159.    Lesen des Datenfiles noetig ist. Dies wird dann aber eher unuebersicht-
  160.    lich.
  161.    Impl.-Hinweis: Reset (filvar, 'filename')
  162.                   oeffnet die Datei 'filename' zum lesen.                    *)
  163.  
  164. PROCEDURE einlesen;
  165.  
  166. VAR
  167.   name    : STRING [25];
  168.   tex     : TEXT;
  169.   n, i, j,
  170.   numcount: INTEGER;
  171.   buff    : STRING [240];
  172.   c1, c2  : CHAR;
  173.  
  174. BEGIN
  175.   pagen;
  176.   Write ('Name des Datenfiles: ');
  177.   ReadLn (name);
  178.   WriteLn;
  179.   Write ('Lese Daten...');
  180.   Reset (tex, name);
  181.   REPEAT                                                     (* BEGIN suchen *)
  182.     ReadLn (tex, buff);
  183.   UNTIL Pos ('BEGIN', buff) = 1;
  184.   spalt := 0;                  (* Spaltenzahl + 1 aus der 1. Zeile bestimmen *)
  185.   c1 := ' ';
  186.   ReadLn (tex, buff);
  187.   FOR i := 1 TO Length (buff) DO
  188.   BEGIN
  189.     c2 := buff[i];
  190.     IF ((c1=' ') AND (c2<>' ')) THEN
  191.       spalt := Succ (spalt);
  192.     c1 := c2;
  193.   END;
  194.   numcount := 0;
  195.   REPEAT
  196.     ReadLn (tex, buff);
  197.     numcount := Succ (numcount);                    (* dabei  Zeilen zaehlen *)
  198.   UNTIL Pos ('END', buff) = 1;
  199.   Reset (tex, name);                                (* Zurueck an den Anfang *)
  200.   REPEAT
  201.     ReadLn (tex, buff);
  202.   UNTIL Pos ('BEGIN', buff) = 1;
  203.   i := 0; n := 0;
  204.   REPEAT                                           (* Koeffizienten einlesen *)
  205.     n := Succ(n);
  206.     FOR i := 1 TO spalt DO
  207.       Read (tex, a^[n,i]);
  208.     Write ('.');
  209.   UNTIL n = numcount;
  210.   zeil := numcount;                                  (* Zeilenzahl festlegen *)
  211.   spalt := Pred (spalt);                          (* Spaltenzahl korrigieren *)
  212.   WriteLn; WriteLn;
  213. END; (* einlesen *)
  214.  
  215. (* ------------------------------------------------------------------------- *)
  216.  
  217. PROCEDURE ende;
  218.  
  219. VAR
  220.   anything : CHAR;
  221.  
  222. BEGIN
  223.   Write ('Mit <RETURN> zurueck zum Desktop');
  224.   Read (anything);
  225. END; (* ende *)
  226.  
  227. (* ------------------------------------------------------------------------- *)
  228.  
  229. BEGIN (* lin_gleich *)
  230.   New (a);
  231.   einlesen;
  232.   householder;
  233.   IF NOT(singmat) THEN
  234.   BEGIN
  235.     ruecktransform;
  236.     ausgabe;
  237.   END;
  238.   ende;
  239. end.
  240.  
  241.  
  242.  
  243.  
  244.  
  245.  
  246.  
  247.