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