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"); }
Design by comdes