/* @(#)glsp.c 17.1.1.1 (ESO-IPG) 01/25/02 17:51:59 */ /* Based on the "Aachen library" glsp.c routines for smooothing splines */ /* ------------------------- MODUL GLSPNP ------------------------- */ /* system includes */ #include #include /* FEROS specific includes */ #include int glspnp #ifdef __STDC__ ( int n, double *xn, double *fn, double *w, int marg_cond, double marg_0, double marg_n, double *a, double *b, double *c, double *d ) #else ( n, xn, fn, w, marg_cond, marg_0, marg_n, a, b, c, d ) int n, marg_cond; double *xn, *fn, *w, marg_0, marg_n, *a, *b, *c, *d; #endif /*********************************************************************** * * * g l s p n p berechnet die Koeffizienten eines nichtparametrischen * * kubischen Ausgleichssplines. Die Art der Randableitung wird durch * * den Paramter marg_cond festgelegt. * * * * Die Splinefunktion wird dargestellt in der Form: * * * * S := S(X) = A(i) + B(i)*(X-XN(i)) + C(i)*(X-XN(i))**2 * * + D(i)*(X-XN(i))**3 * * * * fuer X Element aus [ XN(i) , XN(i+1) ] , i = 0 (1) n-1 . * * * * ==================================================================== * * * * EINGABEPARAMETER: * * ----------------- * * * * Name Typ/Laenge Bedeutung * * ------------------------------------------------------------------- * * n int/-- Index des letzten Knotens * * es gilt: n > 4 fuer marg_cond = 1,2,3 * * n > 5 fuer marg_cond = 4 * * xn double/[n+1] Stuetzstellen des Splines; sie muessen * * str. monoton steigend angeordnet sein * * fn double/[n+1] Messwerte an den Stuetzstellen; * * fuer marg_cond = 4 gilt zusaetzlich: * * fn [0] = fn [n] * * w double/[n+1] Gewichte zu den Messwerten; * * die Gewichte muessen positiv sein; * * fuer marg_cond = 4 gilt zusaetzlich: * * w [0] = w [n] * * marg_cond int/-- Art der Randbedingung * * = 1 : 1. Randableitung vorgegeben * * = 2 : 2. Randableitung vorgegeben * * = 3 : 3. Randableitung vorgegeben * * = 4 : periodische Splinefunktion * * marg_0 double/-- } Randableitung in xn[0] bzw. x[n] * * marg_n double/-- } (nicht von Bedeutung fuer * * marg_cond = 4) * * * * * * AUSGABEPARAMETER: * * ----------------- * * * * Name Typ/Laenge Bedeutung * * ------------------------------------------------------------------- * * a double/[n+1] } Koeffizienten der Splinefunktion * * b double/[n+1] } fuer die Elemente 0 bis n-1 * * c double/[n+1] } ( das n-te Element dient jeweils * * d double/[n+1] } als Hilfselement ) * * * * * * WERT DES UNTERPROGRAMMS: * * ------------------------ * * * * = 0 : kein Fehler * * = 1 : Fehler: Abbruch bei der Zerlegung in den Unterprogrammen * * fdiasy, fdiag, fzyfsy * * = 2 : Fehler: n < 5 bzw. n < 6 * * = 3 : Fehler: die Stuetzstellen sind nicht streng monotom * * = 4 : Fehler: fn [0] ungleich fn [n] oder w [0] ungleich w [n] * * = 5 : Fehler: unzulaessiges Gewicht * * = 6 : Fehler: unzulaessiger Wert fuer marg_cond * * = 7 : Fehler: nicht genuegend Speicherplatz fuer die Hilfsfelder * * * * ==================================================================== * * * * benutzte Unterprogramme: glsp1a, glsp2a, glsp3a, glsppe * * ------------------------ * * * * aus der C-Bibliothek benutzte Unterprogramme: malloc, free * * --------------------------------------------- * * * * benutzte Konstanten: NULL * * -------------------- * * * * ==================================================================== * * * * Bemerkung : Einen natuerlichen Ausgleichsspline erhaelt man mit * * --------- marg_cond = 2 und marg_0 = marg_n = 0.0 * * * ***********************************************************************/ {int error=0, i; double *h, *h1, *h2, *md, *ud1, *ud2, *rs, *hup; /* Plausibilitaetspruefung */ if (n < 5) return (2); for (i=0; i<=n-1; ++i) if (xn [i] >= xn [i+1]) return (3); for (i=0; i<=n; ++i) if (w [i] <= 0.0) return (5); /* Bereitstellung der Hilfsfelder fuer die Unterprogramme */ switch (marg_cond) {case 1: case 2: case 3: {h = (double*) malloc (n * sizeof (double)); if (h == NULL) return (7); h1 = (double*) malloc (n * sizeof (double)); if (h1 == NULL) return (7); h2 = (double*) malloc (n * sizeof (double)); if (h2 == NULL) return (7); md = (double*) malloc (n * sizeof (double)); if (md == NULL) return (7); ud1 = (double*) malloc (n * sizeof (double)); if (ud1 == NULL) return (7); ud2 = (double*) malloc (n * sizeof (double)); if (ud2 == NULL) return (7); rs = (double*) malloc (n * sizeof (double)); if (rs == NULL) return (7); break; } case 4: {h = (double*) malloc ((n+1) * sizeof (double)); if (h == NULL) return (7); h1 = (double*) malloc ((n+1) * sizeof (double)); if (h1 == NULL) return (7); h2 = (double*) malloc ((n+1) * sizeof (double)); if (h2 == NULL) return (7); md = (double*) malloc ((n+1) * sizeof (double)); if (md == NULL) return (7); rs = (double*) malloc ((n+1) * sizeof (double)); if (rs == NULL) return (7); hup = (double*) malloc ((9*n-13) * sizeof (double)); if (hup == NULL) return (7); break; } default: return (6); } /* Aufruf der Unterprogramme, je nach Art der Randbedingung */ /* switch (marg_cond) {case 1: {error = glsp1a (n, xn, fn, w, marg_0, marg_n, 0, a, b, c, d, h, h1, h2, md, ud1, ud2, rs); break;} case 2: {error = glsp2a (n, xn, fn, w, marg_0, marg_n, 0, a, b, c, d, h, h1, h2, md, ud1, ud2, rs); break;} case 3: {error = glsp3a (n, xn, fn, w, marg_0, marg_n, a, b, c, d, h2, md, ud1, ud2, rs, h, h1); break;} case 4: {if (n < 6) return (-1); error = glsppe (n, xn, fn, w, 0, a, b, c, d, h, h1, h2, md, rs, hup); break;} } */ error = glsp2a (n, xn, fn, w, marg_0, marg_n, 0, a, b, c, d, h, h1, h2, md, ud1, ud2, rs); return (error); } /* ------------------------- ENDE GLSPNP ------------------------ */ /* ------------------------- MODUL GLSP2A ------------------------- */ int glsp2a #ifdef __STDC__ ( int n, double *xn, double *fn, double *w, double marg_0, double marg_n, int rep, double *a, double *b, double *c, double *d, double *h, double *h1, double *h2, double *md, double *ud1, double *ud2, double *rs ) #else ( n, xn, fn, w, marg_0, marg_n, rep, a, b, c, d, h, h1, h2, md, ud1, ud2, rs ) int n, rep; double *xn, *fn, *w, marg_0, marg_n, *a, *b, *c, *d; double *h, *h1, *h2, *md, *ud1, *ud2, *rs; #endif /*********************************************************************** * * * g l s p 2 a berechnet die Koeffizienten eines kubischen Aus- * * gleichssplines mit vorgegebener zweiter Randableitung. * * * * Die Splinefunktion wird dargestellt in der Form: * * * * S := S(X) = A(i) + B(i)*(X-XN(i)) + C(i)*(X-XN(i))**2 * * + D(i)*(X-XN(i))**3 * * * * fuer X Element aus [ XN(i) , XN(i+1) ] , i = 0 (1) n-1 . * * * * ==================================================================== * * * * EINGABEPARAMETER: * * ----------------- * * * * Name Typ/Laenge Bedeutung * * ------------------------------------------------------------------- * * n int/-- Index des letzten Knotens * * es gilt: n > 4 * * xn double/[n+1] Stuetzstellen des Splines; sie muessen * * str. monoton steigend angeordnet sein * * fn double/[n+1] Messwerte an den Stuetzstellen * * w double/[n+1] Gewichte zu den Messwerten; * * die Gewichte muessen positiv sein * * marg_0 double/-- 2. Randableitung bei xn[0] * * marg_n double/-- 2. Randableitung bei xn[n] * * rep int/-- gibt an, ob es sich um einen Wieder- * * holungsaufruf handelt: * * = 0 : es handelt sich um den ersten * * Aufruf, d.h. die Matrix zur Be- * * rechnung der c(i) wird aufge- * * baut und mit Hilfe des Unter- * * programms fdiasy zerlegt. * * = 1 : es handelt sich um einen Wie- * * derholungsaufruf, d.h. es wird * * nur die rechte Seite des Glei- * * chungssystems aufgebaut. * * Daher duerfen die Felder md, * * ud1, ud2, h, h1, h2 nach dem * * ersten Aufruf nicht veraendert * * werden! * * * * * * AUSGABEPARAMETER: * * ----------------- * * * * Name Typ/Laenge Bedeutung * * ------------------------------------------------------------------- * * a double/[n+1] } Koeffizienten der Splinefunktion * * b double/[n+1] } fuer die Elemente 0 bis n-1 * * c double/[n+1] } ( Das jeweils n-te Element dient * * d double/[n+1] } als Hilfselement ) * * * * * * HILFSFELDER * * ----------- * * * * Name Typ/Laenge Bedeutung * * ------------------------------------------------------------------- * * h double/[n] } * * h1 double/[n] } Hilfsfelder * * h2 double/[n] } * * md double/[n] } (die ersten Elemente der Felder * * ud1 double/[n] } md, ud1, ud2, rs werden nicht * * ud2 double/[n] } benoetigt) * * rs double/[n] } * * * * * * WERT DES UNTERPROGRAMMS: * * ------------------------ * * * * = 0 : kein Fehler * * = 1 : Fehler: Abbruch in f d i a s y * * = 2 : Fehler: n < 5 * * = 3 : Fehler: falscher Wert fuer rep * * * * ==================================================================== * * * * benutzte Unterprogramme: fdiasy, fdiasl * * ------------------------ * * * * ==================================================================== * * * * Bemerkung : (1) g l s p 2 a sollte von einem Programm aufgerufen * * --------- werden, das die Voraussetzungen ueberprueft (z.B. * * glspnp, glsppa). * * (2) Fuer parametrische Splines mit unterschiedlichen * * Gewichten muss rep auf 0 gesetzt werden. * * (3) Fuer einen natuerlichen Ausgleichsspline waehlt * * man marg_0 = marg_n = 0.0 . * * * ***********************************************************************/ {int i, k, error; double h_var_1, h_var_2; void fdiasl (); if (rep != 0 && rep != 1) return (3); if (!rep) /* Bestimmung der Hilfsgroessen und der Element der Ober- und Hauptdiagonalen des Gleichungssystems */ {for (i=0; i<=n-1; ++i) {h [i] = xn [i+1] - xn [i]; h1 [i] = 1. / h [i]; c [i] = h1 [i] * h1 [i]; b [i] = 6. / w [i]; } b [n] = 6. / w [n]; for (i=0; i<=n-2; ++i) h2 [i] = h1 [i] + h1 [i+1]; /* aeussere Diagonale */ for (i=1; i<=n-3; ++i) ud2 [i] = b [i+1] * h1 [i] * h1 [i+1]; /* innere Diagonale */ for (i=1; i<=n-2; ++i) ud1 [i] = h [i] - b [i] * h1 [i] * h2 [i-1] - b [i+1] * h1 [i] * h2 [i]; /* Hauptdiagonale */ for (i=1; i<=n-1; ++i) {k = i - 1; md [i] = 2. * (h [k] + h [i]) + b [k] * c [k] + b [i] * h2 [k] * h2 [k] + b [i+1] * c [i]; } } /* Berechnung der rechten Seite */ c [0] = 0.5 * marg_0; c [n] = 0.5 * marg_n; h_var_2 = (fn[2] - fn[1]) * h1[1]; h_var_1 = (fn[3] - fn[2]) * h1[2]; rs [1] = 3. * (h_var_2 - (fn[1] - fn[0]) * h1[0]) - c[0] * (h[0] - 6. / w[0] * h1[0] * h1[0] - 6. / w[1] * h1[0] * h2[0]); rs [2] = 3. * (h_var_1 - h_var_2) - c[0] * (6. / w[1]) * h1[0] * h1[1]; for (i=3; i<=n-3; ++i, h_var_1=h_var_2) {h_var_2 = (fn[i+1] - fn[i]) * h1[i]; rs [i] = 3. * (h_var_2 - h_var_1); } h_var_2 = (fn[n-1] - fn[n-2]) * h1[n-2]; rs [n-2] = 3. * (h_var_2 - h_var_1) - c[n] * 6. / w[n-1] * h1[n-2] * h1[n-1]; rs [n-1] = 3. * ((fn[n] - fn[n-1]) * h1[n-1] - h_var_2) - c[n] * (h[n-1] - 6. / w[n-1] * h1[n-1] * h2[n-2] - 6. / w[n] * c[n]); /* Berechnen der Koeffizienten c[i], i=1(1)n-1 durch Loesen des Gleichungssystems mit Hilfe des Unterprogramms FDIASY bzw. FDIASL */ if (! rep) {error = fdiasy (n-1, md, ud1, ud2, rs, c); /* mit vorheriger */ if (error != 0) /* Zerlegung */ if (error == -2) return (2); else return (1); } else /* im Wiederholungsfall ist */ fdiasl (n-1, md, ud1, ud2, rs, c); /* keine Zerlegung notwendig. */ /* Berechnung der restlichen Splinekoeffizienten */ a [0] = fn[0] + 2. / w[0] * h1[0] * (c[0] - c[1]); for (i=1; i<=n-1; ++i) {k = i - 1; a [i] = fn[i] - 2. / w[i] * (c[k] * h1[k] - h2[k] * c[i] + c[i+1] * h1[i]); } a [n] = fn [n] - 2. / w[n] * h1[n-1] * (c[n-1] - c[n]); for (i=0; i<=n-1; ++i) {k = i + 1; b [i] = h1[i] * (a[k] - a[i]) - h[i] / 3. * (c[k] + 2. * c[i]); d [i] = h1[i] / 3. * (c[k] - c[i]); } return (0); } /* ------------------------- ENDE GLSP2A ------------------------ */ /* ------------------------- MODUL FDIAS -------------------------- */ #include int fdiasy #ifdef __STDC__ ( int n, double *md, double *ud1, double *ud2, double *rs, double *x ) #else ( n, md, ud1, ud2, rs, x ) int n; double *md, *ud1, *ud2, *rs, *x; #endif /*********************************************************************** * * * f d i a s y berechnet die L”sung des linearen Gleichungssstems * * A * X = RS mit einer fuenfdiagonalen, symmetrischen, positiv * * definiten Matrix A. * * * * Die Matrix A wird durch die Vektoren md, ud1 und ud2 beschrieben. * * Das Gleichungssystem hat die Form: * * * * md[1]*x[1] + ud1[1]*x[2] + ud2[1]*x[3] = rs[1] * * ud1[1]*x[1] + md[2]*x[2] + ud1[2]*x[3] + ud2[2]*x[4] = rs[2] * * ud2[i-2]*x[i-2] + ud1[i-1]*x[i-1] + md[i]*x[i] * * + ud1[i]*x[i+1] + ud2[i]*x[i+2] = rs[i] * * fuer i=3(1)n-2 * * ud2[n-3]*x[n-2] + ud1[n-2]*x[n-1] + * * + md[n-1]*x[n-1] + ud1[n-1]*x[n] = rs[1] * * ud2[n-2]*x[n-2] + ud1[n-1]*x[n-1] + md[n]*x[n] = rs[1] * * * * * ==================================================================== * * * * EINGABEPARAMETER: * * ----------------- * * * * Name Typ/Laenge Bedeutung * * ------------------------------------------------------------------ * * n int/-- Anzahl der Gleichungen; es gilt: n > 3 * * md double/[n+1] } Elemente der Matrix A * * ud1 double/[n+1] } (das erste Element wird jeweils * * ud2 double/[n+1] } nicht benutzt) * * rs double/[n+1] } md : Hauptdiagonale * * ud1: Diagonale oberhalb der Haupt- * * diagonalen * * ud2: oberste Diagonale * * rs : rechte Seite * * * * * * AUSGABEPARAMETER: * * ----------------- * * * * Name Typ/Laenge Bedeutung * * ------------------------------------------------------------------ * * * * md double/[n+1] } Eingabewerte werden ueberschrieben * * ud1 double/[n+1] } mit Hilfelementen * * ud2 double/[n+1] } * * rs double/[n+1] } * * * * x double/[n+1] Loesung des Gleichungssystems (in den * * Elementen 1 bis n) * * * * * * WERT DES UNTERPROGRAMMS: * * ------------------------ * * * * = 0 : kein Fehler * * = 1 : Fehler: Matrix A ist nicht positiv definit (, da nicht * * streng regulaer) * * = -1: Fehler: Matrix A ist nicht positiv definit * * = -2: Fehler: n < 4 * * * * ==================================================================== * * * * benutzte Unterprogramme: fdiasz, fdiasl * * ------------------------ * * * * ==================================================================== * * * * Bemerkung: Die Determinante von A kann nach dem erfolgreichen * * ---------- Aufruf des Unterprogramms wie folgt berechnet werden: * * det (A) = md[1] * md[2] * ...* md[n] * * * * ==================================================================== * * * * Autorin: Dorothee Seesing-Voelkel * * -------- * * * ***********************************************************************/ {int error; void fdiasl (); if (n < 4) return (-2); error = fdiasz (n, md, ud1, ud2); /* Zerlegung der Matrix */ if (error == 0) fdiasl (n, md, ud1, ud2, rs, x); /* Rueckwaertselimination */ return (error); } #ifdef __STDC__ double sqrt (double x); #else double sqrt (); #endif # #define MACH4_EPS 4.*MACH_EPS int fdiasz #ifdef __STDC__ ( int n, double *md, double *ud1, double *ud2 ) #else ( n, md, ud1, ud2 ) int n; double *md, *ud1, *ud2; #endif /*********************************************************************** * * * f d i a s z zerlegt eine fuenfdiagonale, symmetrische, positiv * * definiten Matrix A in ihre Faktoren B-Transponiert * D * B nach dem * * Cholesky-Verfahren fuer fuenfdiagonale Matrizen. * * (Naehere Beschreibung siehe Modul fdiasy) * * * * ==================================================================== * * * * EINGABEPARAMETER: * * ----------------- * * * * Name Typ/Laenge Bedeutung * * ------------------------------------------------------------------ * * n int/-- Anzahl der Gleichungen; es gilt: n > 3 * * md double/[n+1] } Elemente der Matrix A * * ud1 double/[n+1] } (das erste Element wird jeweils * * ud2 double/[n+1] } nicht benutzt) * * md : Hauptdiagonale * * ud1: Diagonale oberhalb der Haupt- * * diagonalen * * ud2: oberste Diagonale * * * * * * AUSGABEPARAMETER: * * ----------------- * * * * Name Typ/Laenge Bedeutung * * ------------------------------------------------------------------ * * * * md double/[n+1] } Eingabewerte werden ueberschrieben * * ud1 double/[n+1] } mit Feldern, die die Zerlegungs- * * ud2 double/[n+1] } matrix enthalten: * * rs double/[n+1] } ud1, ud2: Diagonalen der normierten * * obere tridiagonale * * Dreiecksmatrix B * * md: Diagonalmatrix D * * * * * * WERT DES UNTERPROGRAMMS: * * ------------------------ * * * * = 0 : kein Fehler * * = 1 : Fehler: Matrix A ist nicht positiv definit (, da nicht * * streng regulaer) * * = -1: Fehler: Matrix A ist nicht positiv definit * * = -2: Fehler: n < 4 * * * * ==================================================================== * * * * benutzte Konstanten: MACH_EPS, MACH4_EPS (Schranke fuer den * * -------------------- relativen Fehler) * * * * ==================================================================== * * * * Autorin: Dorothee Seesing-Voelkel * * -------- * * * ***********************************************************************/ {int error, i; double h_var_1, h_var_2, h_var_3, def_reg; if (n < 4) return (-2); /* Ueberpruefung der Voraussetzung */ /* Fuer n=1 wird die Matrix auf positive Definitheit und strenge Regularitaet ueberprueft */ ud1 [n] = ud2 [n] = ud2 [n-1] = 0.0; def_reg = abs (md[1]) + abs (ud1[1]) + abs (ud2[1]); if (def_reg == 0.0) return (1); def_reg = 1. / def_reg; if (md[1] < 0.0) return (-1); if (abs (md[1]) * def_reg <= MACH4_EPS) return (1); /* Zerlegung der Matrix bei gleichzeitiger Ueberpruefung der positiven Definitheit und der strengen Regularitaet. */ h_var_1 = ud1[1]; ud1 [1] /= md [1]; h_var_2 = ud2[1]; ud2 [1] /= md[1]; def_reg = abs (h_var_1) + abs (md[2]) + abs (ud1[2]) + abs (ud2[2]); if (def_reg == 0.0) return (1); def_reg = 1. / def_reg; md [2] -= h_var_1 * ud1[1]; if (md[2] < 0.0) return (-1); if (abs (md[2]) <= MACH4_EPS) return (1); h_var_1 = ud1[2]; ud1 [2] = (ud1[2] - h_var_2 * ud1[1]) / md[2]; h_var_3 = ud2[2]; ud2 [2] /= md[2]; for (i=3; i<=n; ++i) {def_reg = abs (h_var_2) + abs (h_var_1) + abs (md[i]) + abs (ud1[i]) + abs (ud2[i]); if (def_reg == 0.0) return (1); def_reg = 1. / def_reg; md [i] -= (md[i-1] * sqr(ud1[i-1]) + h_var_2 * ud2[i-2]); if (md[i] < 0.0) return (-1); if (abs (md[i]*def_reg) <= MACH4_EPS) return (1); if (i < n) {h_var_1 = ud1[i]; ud1 [i] = (ud1[i] - h_var_3 * ud1[i-1]) / md[i]; } if (i < n-1) {h_var_2 = h_var_3; h_var_3 = ud2[i]; ud2 [i] /= md[i]; } } return (0); } void fdiasl #ifdef __STDC__ ( int n, double *md, double *ud1, double *ud2, double *rs, double *x ) #else ( n, md, ud1, ud2, rs, x ) int n; double *md, *ud1, *ud2, *rs, *x; #endif /*********************************************************************** * * * f d i a s l berechnet die L”sung des linearen Gleichungssstems * * A * X = RS mit einer fuenfdiagonalen, symmetrischen, positiv * * definiten Matrix A, die in zerlegter Form vorliegt. * * (Naehere Beschreibung siehe Modul fdiasy und fdiasz) * * * * ==================================================================== * * * * EINGABEPARAMETER: * * ----------------- * * * * Name Typ/Laenge Bedeutung * * ------------------------------------------------------------------ * * n int/-- Anzahl der Gleichungen; es gilt: n > 3 * * md double/[n+1] } Zerlegungsmatrix der Matrix A * * ud1 double/[n+1] } (siehe Modul fdiasz) * * ud2 double/[n+1] } md : Diagonalmatrix D * * ud1: } obere Dreiecksmatrix B * * ud2: } * * rs double/[n+1] rs : rechte Seite * * * * * * AUSGABEPARAMETER: * * ----------------- * * * * Name Typ/Laenge Bedeutung * * ------------------------------------------------------------------ * * * * x double/[n+1] Loesung des Gleichungssystems (in den * * Elementen 1 bis n) * * * * ==================================================================== * * * * Autorin: Dorothee Seesing-Voelkel * * -------- * * * ***********************************************************************/ {int i; double h_var_1, h_var_2, h_var_3; /* Vorwaertselimination */ h_var_1 = rs[1]; rs [1] /= md[1]; h_var_2 = rs[2] - ud1[1] * h_var_1; rs [2] = h_var_2 / md[2]; for (i=3; i<=n; ++i) {h_var_1 = rs[i] - ud1[i-1] * h_var_2 - ud2[i-2] * h_var_1; rs [i] = h_var_1 / md[i]; h_var_3 = h_var_2; h_var_2 = h_var_1; h_var_1 = h_var_3; } /* Rueckwaertselimination */ x [n] = rs[n]; x[n-1] = rs[n-1] - ud1[n-1] * x[n]; for (i=n-2; i>=1; --i) x [i] = rs[i] - ud1[i] * x[i+1] - ud2[i] * x[i+2]; return; } /* ------------------------- ENDE FDIAS ------------------------- */ /* ------------------------- MODUL SPVAL -------------------------- */ double spval #ifdef __STDC__ ( int n, double x0, double *a, double *b, double *c, double *d, double *x, double ausg[] ) #else ( n, x0, a, b, c, d, x, ausg ) int n; double x0, ausg [3]; double *a, *b, *c, *d, *x; #endif /*********************************************************************** * * * s p v a l berechnet Funktionswerte eines kubischen Polynom- * * splines und seiner nichttrivialen Ableitungen * * * * ==================================================================== * * * * EINGABEPARAMETER: * * ----------------- * * * * Name Typ/Laenge Bedeutung * * ------------------------------------------------------------------ * * n int/--- Nummer des letzten Knotens * * x0 double/--- Stelle, an der der Funktionswert gesucht * * wird * * a double/n } * * b double/n } Splinekoeffizienten * * c double/n } * * d double/n } * * x double/n+1 Stuetzstellen * * * * AUSGABEPARAMETER: * * ----------------- * * * * Name Typ/Laenge Bedeutung * * ------------------------------------------------------------------ * * ausg double/3 Ableitungen an der Stelle x0 * * ausg [0] enthaelt die 1., * * ausg [1] die 2. und * * ausg [2] die 3. Ableitung * * * * * * WERT DES UNTERPROGRAMMS: * * ------------------------ * * * * Funktionswert an der Stelle x0 * * * * ==================================================================== * * * * Bemerkung: spval muss im rufenden Programm als double verein- * * ---------- bart werden * * * ***********************************************************************/ {int i, k, m; i = 0, k = n; while (m = (i+k) >> 1, m != i) /* das Intervall, in dem x0 */ {if (x0 < x [m]) k = m; /* liegt, wird gesucht. */ else i = m; } x0 -= x [i]; ausg [0] = (3. * d [i] * x0 + 2. * c [i]) * x0 + b [i]; ausg [1] = 6. * d [i] * x0 + 2. * c [i]; ausg [2] = 6. * d [i]; return (((d [i] * x0 + c [i]) * x0 + b [i]) * x0 + a [i]); } /* ------------------------ ENDE SPVAL --------------------------- */