Der folgende Source liefert die Funktionswerte in f[x] für das Intervall eines Splines zwischen dem 4ten und 5ten Intervall. Die Berechnung erfolgt nach der Akima-Interpolation.

/* Für diesen Source benötigen Sie noch den Gauss-Algorithmus zum Lösen der Matrix */
/* Kubische Splines mit Akima-Interpolation */
/* Minimal 2 Punkte, Maximal 5 Punkte (2 vor Intervall,  1 nachher) */
  struct spline {
    WORD von, bis;            /* Normalerweise 0 bis 4, bei Randpunkten entsprechend anders */
    WORD Px[6],Py[6];         /* Maximal 5 Punktepaare */
    BOOL calculated;          /* Schon berechnet? */
    double dq[5];             /* Privat, dqi zur Berechnug von si */
    double si[2];             /* Privat, Steigung */
    double f[4];              /* Funktion f[4]x^3+f[3]x^2+f[1]x+f[0] */
    };

int fakima(struct spline *sp, int intervall)
{
  WORD i;
  WORD d1,d2;
  double a1,a2;

  double *mat[DIM];  /* Zeile,Spalte, Lösungsvektor (Quadratisch) */


  // Test, ob schon berechnet
  if (sp-> calculated == TRUE)
    return(0);

  // Berechnung der dq
  if (sp->von > 2)
  {
    printf("Fehler: Spline braucht min. 2 Stützpunkte!(von)\n");
    return(-1);
  }

  if (sp->bis < 5)
  {
    printf("Fehler: Spline braucht min. 2 Stützpunkte!(bis)\n");
    return(-1);
  }

  // Alle dq's berechnen, die man so berechnen kann
  for (i=sp->von;i < sp->bis;i++)
  {
    d1 = sp->Py[i+1] - sp->Py[i];
    d2 = sp->Px[i+1] - sp->Px[i];
    if (d2 == 0)
      sp->dq[i] = 0.0;
    else
      sp->dq[i] = (double)(d1)/(double)(d2);
  }

  if (sp->von > 1)  // Randpunkt korrigieren
    sp->dq[1] = 2.0 * sp->dq[2] * sp->dq[3];

  if (sp->von > 0)  // Randpunkt korrigieren
    sp->dq[0] = 2.0 * sp->dq[1] * sp->dq[2];

  if (sp->bis < 4)  // Randpunkt korrigieren
    sp->dq[4] = 2.0 * sp->dq[3] * sp->dq[2];

  if (sp->bis < 5)  // Randpunkt korrigieren
    sp->dq[5] = 2.0 * sp->dq[4] * sp->dq[3];


  // si0 berechnen (Steigungen)
  // Startpunkt ist intervall
  a1 = fabs(sp->dq[intervall+1] - sp->dq[intervall]);
  a2 = fabs(sp->dq[intervall-1] - sp->dq[intervall-2]);
  if ((a1+a2) == 0.0)
    sp->si[0] = (sp->dq[intervall-1]+sp->dq[intervall]) / 2;
  else
    sp->si[0] = (a1 * sp->dq[intervall-1] + a2 * sp->dq[intervall]) / (a1+a2);

  // si1 berechnen (Steigungen)
  // Startpunkt ist intervall
  intervall ++;
  a1 = fabs(sp->dq[intervall+1] - sp->dq[intervall]);
  a2 = fabs(sp->dq[intervall-1] - sp->dq[intervall-2]);
  if ((a1+a2) == 0.0)
    sp->si[1] = (sp->dq[intervall-1]+sp->dq[intervall]) / 2;
  else
    sp->si[1] = (a1 * sp->dq[intervall-1] + a2 * sp->dq[intervall]) / (a1+a2);

  // Für Gleichungssystem Speicher
  for (i=0;i<DIM;i++)
    mat[i] = (double *)calloc(DIM,sizeof(double));

  // Gleichungssystem aufbauen
  // Differenzieren

  intervall --;
   mat[0][0] = ((double)(sp->Px[intervall])) * ((double)(sp->Px[intervall])) * ((double)(sp->Px[intervall]));      // x4 einsetzen
   mat[0][1] = ((double)(sp->Px[intervall])) * ((double)(sp->Px[intervall]));
   mat[0][2] = (double)(sp->Px[intervall]);
   mat[0][3] = 1.0;
   mat[0][4] = (double)(sp->Py[intervall]);

   mat[1][0] = ((double)(sp->Px[intervall+1])) * ((double)(sp->Px[intervall+1])) * ((double)(sp->Px[intervall+1]));      // x4 einsetzen
   mat[1][1] = ((double)(sp->Px[intervall+1])) * ((double)(sp->Px[intervall+1]));
   mat[1][2] = (double)(sp->Px[intervall+1]);
   mat[1][3] = 1.0;
   mat[1][4] = (double)(sp->Py[intervall+1]);

   mat[2][0] = 3.0 * ((double)(sp->Px[intervall])) * ((double)(sp->Px[intervall])) ;      // x4 einsetzen
   mat[2][1] = 2.0 * ((double)(sp->Px[intervall])) ;
   mat[2][2] = 1.0;
   mat[2][3] = 0.0;
   mat[2][4] = sp->si[0];

   mat[3][0] = 3.0 * ((double)(sp->Px[intervall+1])) * ((double)(sp->Px[intervall+1]));
   mat[3][1] = 2.0 * ((double)(sp->Px[intervall+1]));
   mat[3][2] = 1.0;
   mat[3][3] = 0.0;
   mat[3][4] = sp->si[1];


   gauss((double **)mat, sp->f,4);

   return(1);
}

int main(int argc, char *argv[])
{
  int i;


/* Kubische Splines mit Akima-Interpolation */
/* Minimal 2 Punkte, Maximal 5 Punkte (2 vor Intervall,  1 nachher) */
{
  struct spline spline;

  spline.von = 0;
  spline.bis = 5;
  spline.Px[0] = -6;
  spline.Py[0] = 4;
  spline.Px[1] = -4;
  spline.Py[1] = 3;
  spline.Px[2] = -1;
  spline.Py[2] = 0;
  spline.Px[3] = 0;
  spline.Py[3] = 0;
  spline.Px[4] = 2;
  spline.Py[4] = 2;
  spline.Px[5] = 3;
  spline.Py[5] = 1;
  spline.calculated = FALSE;

  fakima(&spline, 2);

  for (i=0;i<4;i++)
  {
    printf("%f ",spline.f[i] );
  }
  printf("\n");

}

© 1995 by Thomas Dorn
Design by comdes