home *** CD-ROM | disk | FTP | other *** search
/ Turbo Toolbox / Turbo_Toolbox.iso / spezial / 15 / graphen / reg / spline.pas < prev   
Encoding:
Pascal/Delphi Source File  |  1989-05-18  |  8.0 KB  |  236 lines
  1. { --------------------------------------------------------------- }
  2. {                   File : SPLINE.PAS  (gehört zu WERTE.PAS)      }
  3. {          Copyright (c) : 1989  Heinz Hagemeyer & TOOLBOX        }
  4. {                Sprache : TURBO PASCAL 4.0  (bzw. 3.0)           }
  5. {              Ein Programm zur Auswertung von Meßwerten          }
  6. {                     mit Hilfe von WERTE.PAS                     }
  7. {             (entstanden aus einer alten CPM Routine)            }
  8. { --------------------------------------------------------------- }
  9.  
  10. { Ist die Interpolation mit Hilfe von kubischen Splines möglich   }
  11. {            wird TRUE zurückgegeben, sonst FALSE                 }
  12.  
  13. FUNCTION  Spline_Interpolation (n         : Index      ;
  14.                                 VAR x,y   : Messwerte  ;
  15.                                 VAR Datei : Text      ) : BOOLEAN ;
  16.  
  17. { a_array : hier werden die Koeffizienten des Linearen Gleichungs- }
  18. {           systems gespeichert.                                   }
  19. { x_array : hier kommen die Lösungen hinein.                       }
  20.  
  21. TYPE a_array  = ARRAY [1 .. Spline_Max, 1 .. Spline_Max1] OF REAL;
  22.      x_array  = ARRAY [0 .. Spline_Max                  ] OF REAL;
  23.  
  24. { ------------------------------------------------------------------------ }
  25. { Löst ein Lineares Gleichungssystem nach dem Gauß'schen Algorithmus. Da   }
  26. { es stets bei diesem Problem eine eindeutige Lösung gibt, kann auf eine   }
  27. { Abfrage daraufhin verzichtet werden.                                     }
  28.  
  29. PROCEDURE gauss (n : Index ; VAR x : x_array; VAR a : a_array);
  30.  
  31. VAR i,j,k,l,m : Index ;
  32.     q         : REAL  ;
  33.  
  34. BEGIN
  35.      FOR i := 1 TO pred (n) DO
  36.      BEGIN
  37.           { PIVOT }
  38.           { Suchen des groessten Elementes }
  39.  
  40.           m := i;
  41.           FOR k := Succ (i) TO n DO
  42.               IF Abs (a[k,i]) > Abs (a[m,i]) THEN m := k;
  43.  
  44.           { und Tauschen }
  45.           { IF m <> i THEN  (* Diese Abfrage ist Geschmaksache *) }
  46.           FOR j := 1 TO Succ (n) DO
  47.           BEGIN
  48.                q := a[i,j]; a[i,j] := a [m,j]; a [m,j] := q;
  49.           END;
  50.  
  51.           { Elimination }
  52.  
  53.           FOR k := Succ (i) TO n DO
  54.           BEGIN
  55.                q := a [k,i] / a [i,i] ;
  56.                FOR j := Succ (i) TO  Succ (n) DO
  57.                    a [k,j] := a [k,j] - a [i,j] * q;
  58.           END;
  59.       END;
  60.  
  61.       { Berechnen der Loesungen }
  62.  
  63.       FOR i := n DOWNTO 1 DO
  64.       BEGIN
  65.           q := a[i,succ(n)];
  66.           FOR j := Succ (i) TO n DO q := q - a [i,j] * x [j];
  67.           x [i] := q / a [i,i]
  68.       END;
  69. END;  { Gauss }
  70.  
  71. { -------------------------------------------------------------------------- }
  72. { Mit dieser Procedure wird die Matrix zur Lösung des Lin. Gleichungssystems }
  73. { aus den Meßwerten erzeugt.                                                 }
  74.  
  75. PROCEDURE erzeuge_a (n : Index ; VAR a : a_array; VAR x,y : Messwerte);
  76.  
  77. VAR hi, hi_1 : REAL;
  78.     i,j      : INTEGER;
  79.  
  80. BEGIN
  81.     FOR i := 1 TO n - 2 DO                    { Anfangswerte setzen         }
  82.     FOR j := 1 TO n - 1 DO a [i,j] := 0;
  83.  
  84.     i    := 1;                                { Berechnung der a [i,k] nach }
  85.     hi   := x [i + 2] - x [i + 1 ];           { dem im Begleittext angegebe-}
  86.     hi_1 := x [i + 1] - x [i     ];           { nen Algorithmus             }
  87.  
  88.     a [i,i]     := 2 * (hi_1 + hi);
  89.     a [i,i + 1] := hi;
  90.     a [i,n - 1] := 3 * (   (y [i + 2] - y [i + 1])/hi
  91.                          - (y [i + 1] - y [i    ])/hi_1 );
  92.  
  93.     FOR i := 2 TO n - 2 DO
  94.     BEGIN
  95.         hi   := x [i + 2] - x [i + 1];
  96.         hi_1 := x [i + 1] - x [i    ];
  97.  
  98.         a [i,i - 1] := hi_1;
  99.         a [i,i]     := 2 * (hi_1 + hi);
  100.         a [i,i + 1] := hi;
  101.         a [i,n - 1] := 3 * (   (y [i + 2] - y [i + 1])/hi
  102.                              - (y [i + 1] - y [i    ])/hi_1 );
  103.  
  104.     END;
  105.  
  106.     i := n - 2;
  107.  
  108.     hi   := x [i + 2] - x [i + 1];
  109.     hi_1 := x [i + 1] - x [i    ];
  110.  
  111.     a [i,i - 1] := hi_1;
  112.     a [i,    i] := 2 * (hi_1 + hi);
  113.     a [i,n - 1] := 3 * (   (y [i + 2] - y [i + 1])/hi
  114.                          - (y [i + 1] - y [i    ])/hi_1 );
  115.  
  116. END;  { Erzeuge }
  117.  
  118. { -------------------------------------------------------------------------- }
  119. { Berechnet die Koeffizienten des i-ten Polynom fi(x) = a + bx + cx^2 + bx^3 }
  120.  
  121. PROCEDURE Berechne_Polynom (i   : Index         ;
  122.                             x,y : messwerte     ;
  123.                             ca  : x_array       ;
  124.                             VAR a,b,c,d : REAL) ;
  125.  
  126. VAR hi      : REAL   ;
  127.  
  128. BEGIN
  129.  
  130.      a    := y [i + 1];
  131.      hi   := x [i+2] - x [i+1];
  132.      c    := ca [i];
  133.      b    := (y [i + 2] - a)/hi - hi*( ca [i + 1] + 2* c ) / 3;
  134.      d    := (ca [i + 1] - c ) / (3*hi);
  135.  
  136. END;
  137.  
  138. { -------------------------------------------------------------------------- }
  139. { Schreibt die Daten x und y in die von WERTE.PAS auszuwertenden Datei und   }
  140. { zwar von xa bis xe.                                                        }
  141.  
  142. PROCEDURE Spline_Zug (xa,xe,x0,a,b,c,d : REAL ;
  143.                       VAR t : Text       ) ;
  144.  
  145. VAR x : REAL;
  146.  
  147. FUNCTION y (x : Real) : Real;   { Berechnet f(x) = a + bx + cx^2 + dx^3     }
  148. BEGIN                           { nach dem HORNER - Schema zur Vermeidung   }
  149.                                 { von langsamen Potenzfunktionen.           }
  150.  
  151.      y    :=  a + x*(b + x*( c + d*x)) ;
  152. END;
  153.  
  154. BEGIN   { Spline_Zug }
  155.     x    := xa ;
  156.  
  157.     WHILE x < xe DO
  158.     BEGIN
  159.          WriteLn (t,x,' ',y (x - x0) );
  160.          x := x + Step                ;
  161.     END;
  162.     WriteLn (t,x,' ',y (x - xa) )     ;  { Letzter Zug, nicht vergessen ! }
  163. END;
  164.  
  165. { -------------------------------------------------------------------------- }
  166. { Sotiert nach der Methode der "aufsteigenden Blasen" die Array's x und y.   }
  167. { Ist die langsamste Methode zur Sortierung eines Array's, dafür aber ein-   }
  168. { fach zu programmieren und zu durchschauen. Schnell genug für diesen Zweck. }
  169. { Überprüft gleichzeitig ob ein x - Wert doppelt ist. Dann kann es sich bei  }
  170. { eingegebenen Daten um keine Funktion handeln. Demnach ist eine Interpola-  }
  171. { lation mittels kubischer Splines nicht möglich.                            }
  172.  
  173. FUNCTION bubblesort (n : Index ; VAR x,y : Messwerte) : BOOLEAN;
  174.  
  175. VAR i      : Index  ;
  176.     Merker : REAL   ;
  177.     Sortiert,
  178.     Gleich : BOOLEAN;
  179.  
  180. BEGIN
  181.      REPEAT
  182.            Sortiert := TRUE;
  183.            FOR i := 2 TO n DO
  184.            IF x [i] < x [i - 1]  THEN
  185.            BEGIN
  186.                 Sortiert  := FALSE;
  187.                 Merker    := x [i];
  188.                 x [i]     := x [i - 1];
  189.                 x [i - 1] := Merker;
  190.                 Merker    := y [i];
  191.                 y [i]     := y [i - 1];
  192.                 y [i - 1] := Merker;
  193.            END;
  194.      UNTIL sortiert;
  195.  
  196.      Gleich := TRUE;
  197.      FOR i := 2 TO n DO Gleich := Gleich AND ( x[i] = x[i - 1] );
  198.      Bubblesort := NOT Gleich;
  199. END; { bubblesort }
  200.  
  201. { ------------------------------------------------------------------------ }
  202.  
  203. VAR cc       : ^x_array ;   { Als Zeiger definiert. Spart Platz im Memory  }
  204.     aa       : ^a_array ;   { Bei MS DOS nicht unbedingt erforderlich. Bei }
  205.     i        : Index    ;   { CP/M mit 64k Hauptspeicher war's sehr nötig !}
  206.     a,b,c,d  : REAL;
  207.  
  208. BEGIN  { Spline_Interpolation }
  209.      IF n > Spline_Max THEN
  210.      BEGIN
  211.           Spline_Interpolation := FALSE;
  212.           EXIT;
  213.      END;
  214.      IF (bubblesort (n,x,y)) AND (n>2) THEN
  215.      BEGIN
  216.           Spline_Interpolation := TRUE;
  217.           New (cc); New (aa);
  218.  
  219.           Erzeuge_a (n,aa^,x,y);
  220.           Gauss (n-2,cc^,aa^);
  221.  
  222.           cc^ [0] := 0;
  223.           cc^ [n - 1] := 0;
  224.  
  225.           FOR i := 0 TO n - 2 DO
  226.           BEGIN
  227.                Berechne_Polynom (i,x,y,cc^,a,b,c,d);
  228.                Spline_Zug (x [i + 1] , x [i + 2], x [i + 1],a,b,c,d,Datei);
  229.           END;
  230.           WriteLn (Datei);
  231.           Dispose (cc);
  232.           Dispose (aa);
  233.     END ELSE  Spline_Interpolation := FALSE ;
  234. END; { Spline_Regression }
  235. { ------------------------------------------------------------------------- }
  236. {                          Ende von SPLINE.PAS                              }