Turbo Toolbox/Turbo_Toolbox.iso/spezial/15/linit/linit.pas

Pascal/Delphi Source File | 1988-12-15 | 9.6 KB | 296 lines

PROGRAM Gauss_Seidel;

{ Löst ein lineares Gleichungssystem iterativ nach dem Verfahren von }
{ Gauß - Seidel.       Copyright (c) 1989 DMV - Verlag               }
{                            und Heinz Hagemeyer                     }

USES Crt; { bei Turbo Pascal 3.0 muß diese Zeile entfernt werden     }

CONST Nmax = 50;   { Max. Anzahl der Gleichungen / Unbekannten       }

TYPE Index = 1 .. Nmax;
     Matrix = ARRAY [1 .. Nmax,1 .. Nmax] OF Real;
     Vektor = ARRAY [1 .. Nmax] OF Real;

{ ------------------------------------------------------------------ }

{ TURBO 3.0 benötigt noch diese Procedure :

FUNCTION Readkey : Char;
VAR c : Char;
BEGIN
     Read (kbd,c) ;
     ReadKey := c ;
END;

 ------------------------------------------------------------------- }

FUNCTION Konvergent (VAR A : Matrix; N : Index) : BOOLEAN;

{         Überprüft Spalten - und Zeilensummenkriterium              }
{                     siehe Begleittext                              }

VAR i,k          : Index   ;
    Zeilensumme,
    Spaltensumme : Real    ;
    K1,K2        : BOOLEAN ;
    c            : Char    ;

BEGIN
     K1 := TRUE; K2 := TRUE;          { Annahme : Konvergent        }
     i := 1;

     REPEAT                           { Für jedes i von 1 .. n      }
           Zeilensumme :=  0 ;        { Vorbesetzen der Variablen   }
           Spaltensumme := 0;

           FOR k := 1 TO N DO         { Berechnung der Summen       }
           BEGIN
                Spaltensumme := Spaltensumme + Abs (A[k,i]);
                Zeilensumme  := Zeilensumme  + Abs (A[i,k]);
           END;

           IF A[i,i] = 0 THEN         { Bestimmt nicht konvergent   }
           BEGIN
                K1 := FALSE; K2 := FALSE;
           END ELSE
                               { Bestimmung der Konvergenzkriterien }
           BEGIN
                Spaltensumme := Spaltensumme            ;
                Zeilensumme  := Zeilensumme             ;
                IF Spaltensumme >= Abs (A[i,i] + A[i,i])
                                THEN K1 := FALSE        ;
                IF Zeilensumme  >= Abs (A[i,i] + A[i,i])
                                THEN K2 := FALSE        ;
           END;
           i := Succ (i);              { nächstes i                }
     UNTIL (i > n) OR NOT (K1 OR K2);
                                       { Bestimmung der Konvergenz }

     IF K1 OR K2 THEN Konvergent := TRUE
     ELSE                                 { Warnung ausgeben       }
     BEGIN
       WriteLn ('W A R N U N G   :   Gleichungssytem warscheinlich divergent');
       Write   ('=============       trotzdem versuchen <J/N> ? ');
       REPEAT c := Upcase (Readkey); UNTIL (c='J') OR (c='N');
       Write (c);
       Konvergent :=  c = 'J';       { Kovergenz per Tastatur ! }
     END;
END { Konvergent };

{ ------------------------------------------------------------------------ }

FUNCTION Fehler (VAR Xalt, Xneu : Vektor ; n : Index ) : Real;

{                  Berechnet den Fehler aller  X - Werte                   }
{     kann auf die speziellen Bedürfnisse des Anwenders angepaßt werden    }
{         hier wird die Abweichung der Norm des Vektors berechnet          }

{ ------------------------------------------------------------------------ }

FUNCTION Norm (VAR X : Vektor ; n : Index) : Real ;

{                  Norm des Vektors X                                      }

VAR No : Real ;
     i : Index;

BEGIN
     No := 0;
     FOR i := 1 TO n DO No := No + Sqr (X[i]);
     Norm := Sqrt (No);
END;

BEGIN
     Fehler := Abs (Norm (Xalt,n) - Norm (Xneu,n));
END;

{ -------------------------------------------------------------------------- }

PROCEDURE Ausgabe (VAR X : Vektor; n : Index);
VAR Format : Byte;
{       Format gibt an, ob ein- oder zweispaltig ausgegeben werden soll      }
BEGIN
     IF n > 30 THEN Format := 20 ELSE Format := 60;
     WriteLn;
     FOR n := 1 TO n DO
     BEGIN
          Write ('X':4,n:2,'=':3,X[n]:11:5,'':Format);
     END;
     WriteLn;
END;

{ ------------------------------------------------------------------------- }

PROCEDURE Iteration (VAR A : Matrix; VAR B : Vektor ; n : Index;
                     VAR X : Vektor );

{ 1 Schritt des Verfahrens nach Gauß - Seidel }

{ ------------------------------------------------------------------------- }

FUNCTION ProduktSumme (VAR A : Matrix; VAR X : Vektor; i,n : Index) : Real;
{            Berechnet die Summe Ai,k * Xk, k = 1, .. ,n                    }
VAR k     : Index ;
    Summe : Real  ;

BEGIN
     Summe := 0;
     FOR k := 1 TO n DO Summe := Summe + A[i,k] * X[k];
     ProduktSumme := Summe;
END;

{ ------------------------------------------------------------------------- }

VAR i : Index;

BEGIN   { Iteration }
     FOR i := 1 TO n DO
         X[i] := (B[i] - ProduktSumme (A,X,i,n) + A[i,i] * X[i]) / A[i,i];
                 {               siehe Begleittext                         }
END     { Iteration } ;

{ ------------------------------------------------------------------------ }

PROCEDURE Berechnung (VAR A : Matrix; VAR B : Vektor;
                      VAR X : Vektor; n : Index; i : Byte);

{ Führt i Iterationen durch }
BEGIN
     FOR i := 1 TO i DO Iteration (A,B,n,X);
END;

{ ----------------------------------------------------------------------- }

PROCEDURE BesetzeX (VAR A : Matrix; VAR B,X : Vektor; n : Index );
{   Vorbesetzen des Lösungsvektors X durch Xi := Aii  und i = 1 .. n      }
VAR i : Index;
BEGIN
     FOR i := 1 TO n DO
         IF A [i,i] <> 0
            THEN x[i] := B[i] / A[i,i]
            ELSE
            BEGIN
                 WriteLn ;
                 WriteLn ('Das Gleichungssystem kann mit diesem Verfahren',
                          ' nicht gelöst werden,');
                 WriteLn ('da A [',i:2,',',i:2,'] = 0 ist!');
                 HALT;
            END;
END;

{ ----------------------------------------------------------------------- }

PROCEDURE Eingabe (VAR A : Matrix; VAR B      : Vektor ;
                   VAR n : Index ; VAR Wieoft : Byte    );

{ Die Eingabe kann sowohl über die Tastatur erfolgen, als auch aus einer }
{         Datei, da Pascal die Tastatur wie eine Datei behandelt.        }
CONST  Tab       = 10;

VAR    Datei     : Text;
       Antwort   : Char;
       DateiName : String [20];
       i,k       : Index;
       Tastatur  : Boolean;
       Offset    : Byte;

BEGIN
     ClrScr;
     WriteLn ('Iterative Berechnung eines linearen Gleichungssystem''s nach',
              ' Gauß - Seidel.');
     WriteLn ('============================================================',
              '===============');
     Write ('Dateiname - für Tastatur "CON" - eingeben : ');
     ReadLn (DateiName);
     WriteLn;

     {              Umwandlung in Großbuchstaben zum Vergleichen          }

     FOR i := 1 TO Length (Dateiname)
             DO DateiName [i] := Upcase (DateiName [i]);
     Tastatur := Dateiname = 'CON';
     IF Tastatur THEN Offset := 1 ELSE Offset := 0;

     Assign (Datei, Dateiname);        { Datei versuchsweise öffnen und  }
     {$I-}                             { auf Erfolg überprüfen.          }
          Reset (Datei);
     {$I+}

     IF (IOResult <> 0) THEN
     BEGIN
         WriteLn ('Datei existiert nicht, Programm wird abgebrochen ! ');
         HALT;
     END;

     Write   ('Anzahl der Gleichungen    : '); ReadLn (Datei,n);

    { Wird nicht von der Tastatur sondern aus einer anderen Datei        }
    { gelesen, so werden die eingelesenen Daten auf dem Bildschirm aus-  }
    { gegeben. Die Koeffizienten aber nur, wenn eine Zeile des LGS in    }
    { eine Bildschirmzeile paßt. (siehe weiter unten)                    }

     IF NOT Tastatur THEN Write (n);
     WriteLn;
     IF n < 7 THEN WriteLn ('Und nun die Koeffizienten :');
     WriteLn; WriteLn;

     FOR i := 1 TO N DO
     BEGIN
          FOR k := 1 TO N DO
          BEGIN
               IF Tastatur THEN GotoXY (Tab*(k-1)+1, WhereY - Offset);
               Read (Datei, A [i,k]);
               IF NOT Tastatur AND (n < 7) THEN Write (A [i,k]:Tab:3);
          END;
          IF Tastatur THEN GotoXY (Tab*n+1,WhereY - Offset);
          ReadLn (Datei, B [i]);
          IF NOT Tastatur AND (n<7) THEN WriteLn (B [i]:Tab:3);
          IF     Tastatur THEN WriteLn;
     END;
     Close (Datei);
     WriteLn;
     Write ('Nach wieviel Schritten wünschen Sie eine Lösungsausgabe . ');
     ReadLn (Wieoft);
END; { Einlesen }

{ ----------------------------------------------------------------------- }

FUNCTION Abbruch : Boolean;
{                Wird nach i Iterationen aufgerufen.                      }
VAR c : Char;
BEGIN
     WriteLn;
     Write ('Noch eine Iteration ? <J/N> : ');
     REPEAT c := Upcase (ReadKey); UNTIL (c='J') OR (c='N');
     WriteLn (c:3);
     Abbruch := c = 'N';
END;

{ =================== H A U P T P R O G R A M M =========================== }

VAR
    X,B,Xalt : Vektor  ;
    A        : Matrix  ;
    n        : Index   ;
    WieOft   : Byte    ;

BEGIN
     ClrScr;
     Eingabe (A,B,n,Wieoft);
     BesetzeX (A,B,X,n);
     Window (1,3,80,25);
     ClrScr;
     IF Konvergent (A,n) THEN
     REPEAT
          GotoXY (1,1);
          WriteLn ('---- bitte warten ----');
          Berechnung (A,B,X,n,Wieoft);
          Xalt := X;
          Iteration  (A,B,n,X);
          ClrScr;
          Ausgabe (X,n);
          WriteLn;
          WriteLn ('Abweichung der Norm von Xi = ',Fehler (X,Xalt,n));
    UNTIL Abbruch;
END.