sw/src/pascal/SPLINE.PS

   1 PROGRAM SPLINEINTERPOLATION(INPUT,OUTPUT);
   2
   3   CONST NMAX=25;
   4
   5   VAR   N,I: INTEGER;
   6         X,Y,Z,
   7         LAMBDA,MY,DELTA,P,Q,U: ARRAY[0..NMAX] OF REAL;
   8         S: REAL;
   9         H,
  10         A,B,C,D: ARRAY[1..NMAX] OF REAL;
  11
  12 BEGIN READ(N);  N := N-1;
  13     WRITELN;
  14     WRITELN("S P L I N E - I N T E R P O L A T I O N");
  15     WRITELN;
  16 (*
  17   EINLESEN DER GEGEBENEN PUNKTE, BERECHNUNG DER
  18   LAENGEN H[I] DER EINZELNEN INTERVALLE.
  19 *)
  20     READ(X[0],Y[0]);
  21     FOR I := 1 TO N DO
  22         BEGIN READ(X[I],Y[I]);
  23             H[I] := X[I] - X[I-1]
  24         END;
  25 (*
  26   BERECHNUNG DER KOEFF. LAMBDA[I], MY[I] SOWIE RECHTEN SEITEN DELTA[I]
  27   DES GLEICHUNGSSYSTEMS ZUR BESTIMMUNG DER ZWEITEN ABLEITUNGEN
  28   IN DEN GEGEBENEN PUNKTEN. (TRIDIAGONALMATRIX!)
  29 *)
  30     LAMBDA[0] := 0;  MY[0] := 0;  DELTA[0] := 0;
  31     FOR I := 1 TO N-1 DO
  32         BEGIN S := H[I] + H[I+1];
  33             LAMBDA[I] := H[I+1]/S;  MY[I] := 1 - LAMBDA[I];
  34             DELTA[I] := ( (Y[I+1]-Y[I])/H[I+1] - (Y[I]-Y[I-1])/H[I] )*6/S
  35         END;
  36     LAMBDA[N] := 0;  MY[N] := 0;  DELTA[N] := 0;
  37 (*
  38   BERECHNUNG DER ZWEITEN ABLEITUNGEN Z[I]
  39 *)
  40     Q[0] := -LAMBDA[0]/2;  U[0] := DELTA[0]/2;
  41     FOR I := 1 TO N DO
  42         BEGIN
  43             P[I] := MY[I]*Q[I-1] + 2;
  44             Q[I] := -LAMBDA[I]/P[I];
  45             U[I] := (DELTA[I] - MY[I]*U[I-1])/P[I]
  46         END;
  47     Z[N] := U[N];
  48     FOR I := N-1 DOWNTO 0 DO Z[I] := Q[I]*Z[I+1] + U[I];
  49 (*
  50   BERECHNUNG DER KOEFFIZIENTEN DER KUBISCHEN POLYNOME
  51         A[I]*X^3 + B[I]*X^2 + C[I]*X + D[I]
  52   FUER DIE EINZELNEN INTERVALLE.
  53 *)
  54     FOR I := 1 TO N DO
  55         BEGIN S := X[I] - X[I-1];
  56             A[I] := (Z[I]-Z[I-1])/(6*S);
  57             B[I] := (Z[I]-6*A[I]*X[I])/2;
  58             C[I] := (Y[I]-Y[I-1])/S -
  59                     A[I]*(SQR(X[I]) + X[I]*X[I-1] + SQR(X[I-1])) -
  60                     B[I]*(X[I]+X[I-1]);
  61             D[I] := Y[I] - ((A[I]*X[I]+B[I])*X[I]+C[I])*X[I]
  62         END;
  63 (*
  64   AUSGABE
  65 *)
  66     WRITELN("P 0:", X[0] :8:2, Y[0] :8:2);
  67     FOR I := 1 TO N DO
  68         BEGIN
  69             WRITELN(" ":5, A[I],B[I],C[I],D[I]);
  70             WRITELN("P", I:2, ":", X[I] :8:2, Y[I] :8:2)
  71         END;
  72     WRITELN
  73 END.