2
/* ------------------------- MODUL GLSPNP ------------------------- */
6
Based on the "Aachen library"
10
routines for smooothing splines
23
/* FEROS specific includes */
31
(int n, double *xn, double *fn, double *w,
32
int marg_cond, double marg_0, double marg_n,
33
double *a, double *b, double *c, double *d)
35
(n, xn, fn, w, marg_cond, marg_0, marg_n, a, b, c, d)
37
double *xn, *fn, *w, marg_0, marg_n, *a, *b, *c, *d;
40
/***********************************************************************
42
* g l s p n p berechnet die Koeffizienten eines nichtparametrischen *
43
* kubischen Ausgleichssplines. Die Art der Randableitung wird durch *
44
* den Paramter marg_cond festgelegt. *
46
* Die Splinefunktion wird dargestellt in der Form: *
48
* S := S(X) = A(i) + B(i)*(X-XN(i)) + C(i)*(X-XN(i))**2 *
49
* + D(i)*(X-XN(i))**3 *
51
* fuer X Element aus [ XN(i) , XN(i+1) ] , i = 0 (1) n-1 . *
53
* ==================================================================== *
58
* Name Typ/Laenge Bedeutung *
59
* ------------------------------------------------------------------- *
60
* n int/-- Index des letzten Knotens *
61
* es gilt: n > 4 fuer marg_cond = 1,2,3 *
62
* n > 5 fuer marg_cond = 4 *
63
* xn double/[n+1] Stuetzstellen des Splines; sie muessen *
64
* str. monoton steigend angeordnet sein *
65
* fn double/[n+1] Messwerte an den Stuetzstellen; *
66
* fuer marg_cond = 4 gilt zusaetzlich: *
68
* w double/[n+1] Gewichte zu den Messwerten; *
69
* die Gewichte muessen positiv sein; *
70
* fuer marg_cond = 4 gilt zusaetzlich: *
72
* marg_cond int/-- Art der Randbedingung *
73
* = 1 : 1. Randableitung vorgegeben *
74
* = 2 : 2. Randableitung vorgegeben *
75
* = 3 : 3. Randableitung vorgegeben *
76
* = 4 : periodische Splinefunktion *
77
* marg_0 double/-- } Randableitung in xn[0] bzw. x[n] *
78
* marg_n double/-- } (nicht von Bedeutung fuer *
85
* Name Typ/Laenge Bedeutung *
86
* ------------------------------------------------------------------- *
87
* a double/[n+1] } Koeffizienten der Splinefunktion *
88
* b double/[n+1] } fuer die Elemente 0 bis n-1 *
89
* c double/[n+1] } ( das n-te Element dient jeweils *
90
* d double/[n+1] } als Hilfselement ) *
93
* WERT DES UNTERPROGRAMMS: *
94
* ------------------------ *
97
* = 1 : Fehler: Abbruch bei der Zerlegung in den Unterprogrammen *
98
* fdiasy, fdiag, fzyfsy *
99
* = 2 : Fehler: n < 5 bzw. n < 6 *
100
* = 3 : Fehler: die Stuetzstellen sind nicht streng monotom *
101
* = 4 : Fehler: fn [0] ungleich fn [n] oder w [0] ungleich w [n] *
102
* = 5 : Fehler: unzulaessiges Gewicht *
103
* = 6 : Fehler: unzulaessiger Wert fuer marg_cond *
104
* = 7 : Fehler: nicht genuegend Speicherplatz fuer die Hilfsfelder *
106
* ==================================================================== *
108
* benutzte Unterprogramme: glsp1a, glsp2a, glsp3a, glsppe *
109
* ------------------------ *
111
* aus der C-Bibliothek benutzte Unterprogramme: malloc, free *
112
* --------------------------------------------- *
114
* benutzte Konstanten: NULL *
115
* -------------------- *
117
* ==================================================================== *
119
* Bemerkung : Einen natuerlichen Ausgleichsspline erhaelt man mit *
120
* --------- marg_cond = 2 und marg_0 = marg_n = 0.0 *
122
***********************************************************************/
126
double *h, *h1, *h2, *md, *ud1, *ud2, *rs, *hup;
129
ud1 = ud2 = (double *) 0;
132
Plausibilitaetspruefung
134
if (n < 5) return (2);
135
for (i = 0; i <= n - 1; ++i)
136
if (xn[i] >= xn[i + 1]) return (3);
137
for (i = 0; i <= n; ++i)
138
if (w[i] <= 0.0) return (5);
140
Bereitstellung der Hilfsfelder fuer die Unterprogramme
148
h = (double *) malloc ((size_t)(n*sizeof (double)));
151
h1 = (double *) malloc ((size_t)(n*sizeof (double)));
154
h2 = (double *) malloc ((size_t)(n*sizeof (double)));
157
md = (double *) malloc ((size_t)(n * sizeof (double)));
160
ud1 = (double *) malloc ((size_t)(n * sizeof (double)));
163
ud2 = (double *) malloc ((size_t)(n * sizeof (double)));
166
rs = (double *) malloc ((size_t)(n * sizeof (double)));
173
h = (double *) malloc ((size_t)((n + 1) * sizeof (double)));
176
h1 = (double *) malloc ((size_t)((n + 1) * sizeof (double)));
179
h2 = (double *) malloc ((size_t)((n + 1) * sizeof (double)));
182
md = (double *) malloc ((size_t)((n + 1) * sizeof (double)));
185
rs = (double *) malloc ((size_t)((n + 1) * sizeof (double)));
188
hup = (double *) malloc ((size_t)((9 * n - 13) * sizeof (double)));
197
Aufruf der Unterprogramme, je nach Art der Randbedingung
202
{case 1: {error = glsp1a (n, xn, fn, w, marg_0, marg_n, 0,
203
a, b, c, d, h, h1, h2, md, ud1, ud2, rs);
205
case 2: {error = glsp2a (n, xn, fn, w, marg_0, marg_n, 0,
206
a, b, c, d, h, h1, h2, md, ud1, ud2, rs);
208
case 3: {error = glsp3a (n, xn, fn, w, marg_0, marg_n, a, b, c, d,
209
h2, md, ud1, ud2, rs, h, h1);
211
case 4: {if (n < 6) return (-1);
212
error = glsppe (n, xn, fn, w, 0, a, b, c, d,
213
h, h1, h2, md, rs, hup);
217
error = glsp2a (n, xn, fn, w, marg_0, marg_n, 0,
218
a, b, c, d, h, h1, h2, md, ud1, ud2, rs);
222
/* ------------------------- ENDE GLSPNP ------------------------ */
223
/* ------------------------- MODUL GLSP2A ------------------------- */
229
(int n, double *xn, double *fn, double *w,
230
double marg_0, double marg_n, int rep,
231
double *a, double *b, double *c, double *d,
232
double *h, double *h1, double *h2, double *md,
233
double *ud1, double *ud2, double *rs)
237
(n, xn, fn, w, marg_0, marg_n, rep, a, b, c, d,
238
h, h1, h2, md, ud1, ud2, rs)
240
double *xn, *fn, *w, marg_0, marg_n, *a, *b, *c, *d;
241
double *h, *h1, *h2, *md, *ud1, *ud2, *rs;
243
/***********************************************************************
245
* g l s p 2 a berechnet die Koeffizienten eines kubischen Aus- *
246
* gleichssplines mit vorgegebener zweiter Randableitung. *
248
* Die Splinefunktion wird dargestellt in der Form: *
250
* S := S(X) = A(i) + B(i)*(X-XN(i)) + C(i)*(X-XN(i))**2 *
251
* + D(i)*(X-XN(i))**3 *
253
* fuer X Element aus [ XN(i) , XN(i+1) ] , i = 0 (1) n-1 . *
255
* ==================================================================== *
257
* EINGABEPARAMETER: *
258
* ----------------- *
260
* Name Typ/Laenge Bedeutung *
261
* ------------------------------------------------------------------- *
262
* n int/-- Index des letzten Knotens *
264
* xn double/[n+1] Stuetzstellen des Splines; sie muessen *
265
* str. monoton steigend angeordnet sein *
266
* fn double/[n+1] Messwerte an den Stuetzstellen *
267
* w double/[n+1] Gewichte zu den Messwerten; *
268
* die Gewichte muessen positiv sein *
269
* marg_0 double/-- 2. Randableitung bei xn[0] *
270
* marg_n double/-- 2. Randableitung bei xn[n] *
271
* rep int/-- gibt an, ob es sich um einen Wieder- *
272
* holungsaufruf handelt: *
273
* = 0 : es handelt sich um den ersten *
274
* Aufruf, d.h. die Matrix zur Be- *
275
* rechnung der c(i) wird aufge- *
276
* baut und mit Hilfe des Unter- *
277
* programms fdiasy zerlegt. *
278
* = 1 : es handelt sich um einen Wie- *
279
* derholungsaufruf, d.h. es wird *
280
* nur die rechte Seite des Glei- *
281
* chungssystems aufgebaut. *
282
* Daher duerfen die Felder md, *
283
* ud1, ud2, h, h1, h2 nach dem *
284
* ersten Aufruf nicht veraendert *
288
* AUSGABEPARAMETER: *
289
* ----------------- *
291
* Name Typ/Laenge Bedeutung *
292
* ------------------------------------------------------------------- *
293
* a double/[n+1] } Koeffizienten der Splinefunktion *
294
* b double/[n+1] } fuer die Elemente 0 bis n-1 *
295
* c double/[n+1] } ( Das jeweils n-te Element dient *
296
* d double/[n+1] } als Hilfselement ) *
302
* Name Typ/Laenge Bedeutung *
303
* ------------------------------------------------------------------- *
305
* h1 double/[n] } Hilfsfelder *
307
* md double/[n] } (die ersten Elemente der Felder *
308
* ud1 double/[n] } md, ud1, ud2, rs werden nicht *
309
* ud2 double/[n] } benoetigt) *
313
* WERT DES UNTERPROGRAMMS: *
314
* ------------------------ *
316
* = 0 : kein Fehler *
317
* = 1 : Fehler: Abbruch in f d i a s y *
318
* = 2 : Fehler: n < 5 *
319
* = 3 : Fehler: falscher Wert fuer rep *
321
* ==================================================================== *
323
* benutzte Unterprogramme: fdiasy, fdiasl *
324
* ------------------------ *
326
* ==================================================================== *
328
* Bemerkung : (1) g l s p 2 a sollte von einem Programm aufgerufen *
329
* --------- werden, das die Voraussetzungen ueberprueft (z.B. *
331
* (2) Fuer parametrische Splines mit unterschiedlichen *
332
* Gewichten muss rep auf 0 gesetzt werden. *
333
* (3) Fuer einen natuerlichen Ausgleichsspline waehlt *
334
* man marg_0 = marg_n = 0.0 . *
336
***********************************************************************/
340
double h_var_1, h_var_2;
343
if (rep != 0 && rep != 1)
346
/* Bestimmung der Hilfsgroessen und der Element der
347
Ober- und Hauptdiagonalen des Gleichungssystems */
349
for (i = 0; i <= n - 1; ++i)
351
h[i] = xn[i + 1] - xn[i];
353
c[i] = h1[i] * h1[i];
358
for (i = 0; i <= n - 2; ++i)
359
h2[i] = h1[i] + h1[i + 1];
363
for (i = 1; i <= n - 3; ++i)
364
ud2[i] = b[i + 1] * h1[i] * h1[i + 1];
368
for (i = 1; i <= n - 2; ++i)
369
ud1[i] = h[i] - b[i] * h1[i] * h2[i - 1] - b[i + 1] * h1[i] * h2[i];
373
for (i = 1; i <= n - 1; ++i)
376
md[i] = 2. * (h[k] + h[i]) + b[k] * c[k]
377
+ b[i] * h2[k] * h2[k] + b[i + 1] * c[i];
381
Berechnung der rechten Seite
386
h_var_2 = (fn[2] - fn[1]) * h1[1];
387
h_var_1 = (fn[3] - fn[2]) * h1[2];
388
rs[1] = 3. * (h_var_2 - (fn[1] - fn[0]) * h1[0]) - c[0] *
389
(h[0] - 6. / w[0] * h1[0] * h1[0] - 6. / w[1] * h1[0] * h2[0]);
390
rs[2] = 3. * (h_var_1 - h_var_2) - c[0] * (6. / w[1]) * h1[0] * h1[1];
391
for (i = 3; i <= n - 3; ++i, h_var_1 = h_var_2)
393
h_var_2 = (fn[i + 1] - fn[i]) * h1[i];
394
rs[i] = 3. * (h_var_2 - h_var_1);
396
h_var_2 = (fn[n - 1] - fn[n - 2]) * h1[n - 2];
397
rs[n - 2] = 3. * (h_var_2 - h_var_1)
398
- c[n] * 6. / w[n - 1] * h1[n - 2] * h1[n - 1];
399
rs[n - 1] = 3. * ((fn[n] - fn[n - 1]) * h1[n - 1] - h_var_2)
400
- c[n] * (h[n - 1] - 6. / w[n - 1] * h1[n - 1] * h2[n - 2]
403
Berechnen der Koeffizienten c[i], i=1(1)n-1
404
durch Loesen des Gleichungssystems mit Hilfe
405
des Unterprogramms FDIASY bzw. FDIASL
409
error = fdiasy (n - 1, md, ud1, ud2, rs, c); /* mit vorheriger */
410
if (error != 0) /* Zerlegung */
418
else /* im Wiederholungsfall ist */
419
fdiasl (n - 1, md, ud1, ud2, rs, c); /* keine Zerlegung notwendig. */
421
Berechnung der restlichen Splinekoeffizienten
423
a[0] = fn[0] + 2. / w[0] * h1[0] * (c[0] - c[1]);
424
for (i = 1; i <= n - 1; ++i)
427
a[i] = fn[i] - 2. / w[i] * (c[k] * h1[k] - h2[k] * c[i]
430
a[n] = fn[n] - 2. / w[n] * h1[n - 1] * (c[n - 1] - c[n]);
432
for (i = 0; i <= n - 1; ++i)
435
b[i] = h1[i] * (a[k] - a[i]) - h[i] / 3. * (c[k] + 2. * c[i]);
436
d[i] = h1[i] / 3. * (c[k] - c[i]);
442
/* ------------------------- ENDE GLSP2A ------------------------ */
443
/* ------------------------- MODUL FDIAS -------------------------- */
448
(int n, double *md, double *ud1, double *ud2, double *rs, double *x)
450
(n, md, ud1, ud2, rs, x)
452
double *md, *ud1, *ud2, *rs, *x;
454
/***********************************************************************
456
* f d i a s y berechnet die L�sung des linearen Gleichungssstems *
457
* A * X = RS mit einer fuenfdiagonalen, symmetrischen, positiv *
458
* definiten Matrix A. *
460
* Die Matrix A wird durch die Vektoren md, ud1 und ud2 beschrieben. *
461
* Das Gleichungssystem hat die Form: *
463
* md[1]*x[1] + ud1[1]*x[2] + ud2[1]*x[3] = rs[1] *
464
* ud1[1]*x[1] + md[2]*x[2] + ud1[2]*x[3] + ud2[2]*x[4] = rs[2] *
465
* ud2[i-2]*x[i-2] + ud1[i-1]*x[i-1] + md[i]*x[i] *
466
* + ud1[i]*x[i+1] + ud2[i]*x[i+2] = rs[i] *
468
* ud2[n-3]*x[n-2] + ud1[n-2]*x[n-1] + *
469
* + md[n-1]*x[n-1] + ud1[n-1]*x[n] = rs[1] *
470
* ud2[n-2]*x[n-2] + ud1[n-1]*x[n-1] + md[n]*x[n] = rs[1] *
472
* ==================================================================== *
474
* EINGABEPARAMETER: *
475
* ----------------- *
477
* Name Typ/Laenge Bedeutung *
478
* ------------------------------------------------------------------ *
479
* n int/-- Anzahl der Gleichungen; es gilt: n > 3 *
480
* md double/[n+1] } Elemente der Matrix A *
481
* ud1 double/[n+1] } (das erste Element wird jeweils *
482
* ud2 double/[n+1] } nicht benutzt) *
483
* rs double/[n+1] } md : Hauptdiagonale *
484
* ud1: Diagonale oberhalb der Haupt- *
486
* ud2: oberste Diagonale *
487
* rs : rechte Seite *
490
* AUSGABEPARAMETER: *
491
* ----------------- *
493
* Name Typ/Laenge Bedeutung *
494
* ------------------------------------------------------------------ *
496
* md double/[n+1] } Eingabewerte werden ueberschrieben *
497
* ud1 double/[n+1] } mit Hilfelementen *
498
* ud2 double/[n+1] } *
499
* rs double/[n+1] } *
501
* x double/[n+1] Loesung des Gleichungssystems (in den *
502
* Elementen 1 bis n) *
505
* WERT DES UNTERPROGRAMMS: *
506
* ------------------------ *
508
* = 0 : kein Fehler *
509
* = 1 : Fehler: Matrix A ist nicht positiv definit (, da nicht *
511
* = -1: Fehler: Matrix A ist nicht positiv definit *
512
* = -2: Fehler: n < 4 *
514
* ==================================================================== *
516
* benutzte Unterprogramme: fdiasz, fdiasl *
517
* ------------------------ *
519
* ==================================================================== *
521
* Bemerkung: Die Determinante von A kann nach dem erfolgreichen *
522
* ---------- Aufruf des Unterprogramms wie folgt berechnet werden: *
523
* det (A) = md[1] * md[2] * ...* md[n] *
525
* ==================================================================== *
527
* Autorin: Dorothee Seesing-Voelkel *
530
***********************************************************************/
539
error = fdiasz (n, md, ud1, ud2); /* Zerlegung der Matrix */
541
fdiasl (n, md, ud1, ud2, rs, x); /* Rueckwaertselimination */
547
double sqrt (double x);
552
#define MACH4_EPS 4.*MACH_EPS
556
(int n, double *md, double *ud1, double *ud2)
560
double *md, *ud1, *ud2;
562
/***********************************************************************
564
* f d i a s z zerlegt eine fuenfdiagonale, symmetrische, positiv *
565
* definiten Matrix A in ihre Faktoren B-Transponiert * D * B nach dem *
566
* Cholesky-Verfahren fuer fuenfdiagonale Matrizen. *
567
* (Naehere Beschreibung siehe Modul fdiasy) *
569
* ==================================================================== *
571
* EINGABEPARAMETER: *
572
* ----------------- *
574
* Name Typ/Laenge Bedeutung *
575
* ------------------------------------------------------------------ *
576
* n int/-- Anzahl der Gleichungen; es gilt: n > 3 *
577
* md double/[n+1] } Elemente der Matrix A *
578
* ud1 double/[n+1] } (das erste Element wird jeweils *
579
* ud2 double/[n+1] } nicht benutzt) *
580
* md : Hauptdiagonale *
581
* ud1: Diagonale oberhalb der Haupt- *
583
* ud2: oberste Diagonale *
586
* AUSGABEPARAMETER: *
587
* ----------------- *
589
* Name Typ/Laenge Bedeutung *
590
* ------------------------------------------------------------------ *
592
* md double/[n+1] } Eingabewerte werden ueberschrieben *
593
* ud1 double/[n+1] } mit Feldern, die die Zerlegungs- *
594
* ud2 double/[n+1] } matrix enthalten: *
595
* rs double/[n+1] } ud1, ud2: Diagonalen der normierten *
596
* obere tridiagonale *
598
* md: Diagonalmatrix D *
601
* WERT DES UNTERPROGRAMMS: *
602
* ------------------------ *
604
* = 0 : kein Fehler *
605
* = 1 : Fehler: Matrix A ist nicht positiv definit (, da nicht *
607
* = -1: Fehler: Matrix A ist nicht positiv definit *
608
* = -2: Fehler: n < 4 *
610
* ==================================================================== *
612
* benutzte Konstanten: MACH_EPS, MACH4_EPS (Schranke fuer den *
613
* -------------------- relativen Fehler) *
615
* ==================================================================== *
617
* Autorin: Dorothee Seesing-Voelkel *
620
***********************************************************************/
624
double h_var_1, h_var_2, h_var_3, def_reg;
627
return (-2); /* Ueberpruefung der Voraussetzung */
629
Fuer n=1 wird die Matrix auf positive Definitheit
630
und strenge Regularitaet ueberprueft
632
ud1[n] = ud2[n] = ud2[n - 1] = 0.0;
633
def_reg = abs (md[1]) + abs (ud1[1]) + abs (ud2[1]);
636
def_reg = 1. / def_reg;
639
if (abs (md[1]) * def_reg <= MACH4_EPS)
642
Zerlegung der Matrix bei gleichzeitiger Ueberpruefung
643
der positiven Definitheit und der strengen Regularitaet.
649
def_reg = abs (h_var_1) + abs (md[2]) + abs (ud1[2]) + abs (ud2[2]);
652
def_reg = 1. / def_reg;
653
md[2] -= h_var_1 * ud1[1];
656
if (abs (md[2]) <= MACH4_EPS)
660
ud1[2] = (ud1[2] - h_var_2 * ud1[1]) / md[2];
663
for (i = 3; i <= n; ++i)
665
def_reg = abs (h_var_2) + abs (h_var_1) + abs (md[i])
666
+ abs (ud1[i]) + abs (ud2[i]);
669
def_reg = 1. / def_reg;
671
md[i] -= (md[i - 1] * sqr (ud1[i - 1]) + h_var_2 * ud2[i - 2]);
674
if (abs (md[i] * def_reg) <= MACH4_EPS)
679
ud1[i] = (ud1[i] - h_var_3 * ud1[i - 1]) / md[i];
693
(int n, double *md, double *ud1, double *ud2, double *rs, double *x)
695
(n, md, ud1, ud2, rs, x)
697
double *md, *ud1, *ud2, *rs, *x;
699
/***********************************************************************
701
* f d i a s l berechnet die L�sung des linearen Gleichungssstems *
702
* A * X = RS mit einer fuenfdiagonalen, symmetrischen, positiv *
703
* definiten Matrix A, die in zerlegter Form vorliegt. *
704
* (Naehere Beschreibung siehe Modul fdiasy und fdiasz) *
706
* ==================================================================== *
708
* EINGABEPARAMETER: *
709
* ----------------- *
711
* Name Typ/Laenge Bedeutung *
712
* ------------------------------------------------------------------ *
713
* n int/-- Anzahl der Gleichungen; es gilt: n > 3 *
714
* md double/[n+1] } Zerlegungsmatrix der Matrix A *
715
* ud1 double/[n+1] } (siehe Modul fdiasz) *
716
* ud2 double/[n+1] } md : Diagonalmatrix D *
717
* ud1: } obere Dreiecksmatrix B *
719
* rs double/[n+1] rs : rechte Seite *
722
* AUSGABEPARAMETER: *
723
* ----------------- *
725
* Name Typ/Laenge Bedeutung *
726
* ------------------------------------------------------------------ *
728
* x double/[n+1] Loesung des Gleichungssystems (in den *
729
* Elementen 1 bis n) *
731
* ==================================================================== *
733
* Autorin: Dorothee Seesing-Voelkel *
736
***********************************************************************/
740
double h_var_1, h_var_2, h_var_3;
746
h_var_2 = rs[2] - ud1[1] * h_var_1;
747
rs[2] = h_var_2 / md[2];
748
for (i = 3; i <= n; ++i)
750
h_var_1 = rs[i] - ud1[i - 1] * h_var_2 - ud2[i - 2] * h_var_1;
751
rs[i] = h_var_1 / md[i];
757
Rueckwaertselimination
760
x[n - 1] = rs[n - 1] - ud1[n - 1] * x[n];
761
for (i = n - 2; i >= 1; --i)
762
x[i] = rs[i] - ud1[i] * x[i + 1] - ud2[i] * x[i + 2];
767
/* ------------------------- ENDE FDIAS ------------------------- */
768
/* ------------------------- MODUL SPVAL -------------------------- */
774
(int n, double x0, double *a, double *b, double *c,
775
double *d, double *x, double ausg[])
777
(n, x0, a, b, c, d, x, ausg)
780
double *a, *b, *c, *d, *x;
782
/***********************************************************************
784
* s p v a l berechnet Funktionswerte eines kubischen Polynom- *
785
* splines und seiner nichttrivialen Ableitungen *
787
* ==================================================================== *
789
* EINGABEPARAMETER: *
790
* ----------------- *
792
* Name Typ/Laenge Bedeutung *
793
* ------------------------------------------------------------------ *
794
* n int/--- Nummer des letzten Knotens *
795
* x0 double/--- Stelle, an der der Funktionswert gesucht *
798
* b double/n } Splinekoeffizienten *
801
* x double/n+1 Stuetzstellen *
803
* AUSGABEPARAMETER: *
804
* ----------------- *
806
* Name Typ/Laenge Bedeutung *
807
* ------------------------------------------------------------------ *
808
* ausg double/3 Ableitungen an der Stelle x0 *
809
* ausg [0] enthaelt die 1., *
810
* ausg [1] die 2. und *
811
* ausg [2] die 3. Ableitung *
814
* WERT DES UNTERPROGRAMMS: *
815
* ------------------------ *
817
* Funktionswert an der Stelle x0 *
819
* ==================================================================== *
821
* Bemerkung: spval muss im rufenden Programm als double verein- *
822
* ---------- bart werden *
824
***********************************************************************/
830
while (m = (i + k) >> 1, m != i) /* das Intervall, in dem x0 */
833
k = m; /* liegt, wird gesucht. */
838
ausg[0] = (3. * d[i] * x0 + 2. * c[i]) * x0 + b[i];
839
ausg[1] = 6. * d[i] * x0 + 2. * c[i];
842
return (((d[i] * x0 + c[i]) * x0 + b[i]) * x0 + a[i]);
845
/* ------------------------ ENDE SPVAL --------------------------- */
added
removed
stdred/feros/libsrc/glsp.c
