periodic cspline first derivative discontinuity
David Necas (Yeti)
yeti@physics.muni.cz
Sat Apr 6 04:05:00 GMT 2002
Hello,
I've found a bug (or I think I've found a bug ;)
in periodic cubic spline implementation  its
first derivative is not continuous at the edges of
the first segment.
This is because length of next segment (h_ip1
in the patch below) is computed modulo even for
the last segment and so it results as negative
(well, not much clearly explained, perhaps the
patch itself is clearer).
Please note applying this patch actually fixes
only first derivative discontinuity in right edge
of the first segment, some dicontinuity may still
appear at the left edge due to probems in cyclic
tridiagonal linear system solver  see my next
bugreport.
Yeti
 gsl1.1.1.orig/interpolation/cspline.c Sun Mar 31 19:32:20 2002
+++ gsl1.1.1/interpolation/cspline.c Mon Apr 1 14:52:12 2002
@@ 168,13 +168,25 @@
return GSL_SUCCESS;
} else {
 state>offdiag[max_index] = xa[1]xa[0];

 for (i = 0; i < sys_size; i++) {
 const double h_i = xa[i + 1]  xa[i];
 const double h_ip1 = xa[(i + 2) % num_points]  xa[i + 1];
+ for (i = 0; i < sys_size1; i++) {
+ const double h_i = xa[i + 1]  xa[i];
+ const double h_ip1 = xa[i + 2]  xa[i + 1];
+ const double ydiff_i = ya[i + 1]  ya[i];
+ const double ydiff_ip1 = ya[i + 2]  ya[i + 1];
+ state>offdiag[i] = h_ip1;
+ state>diag[i] = 2.0 * (h_ip1 + h_i);
+ state>g[i] = 3.0 * (ydiff_ip1 / h_ip1  ydiff_i / h_i);
+ }
+ /* fixed derivative discontinuity in xa[1], but not in xa[0] 
+ * this one is due to some tridiag solver problem
+ * FIXME
+ */
+ {
+ i = sys_size  1;
+ const double h_i = xa[i + 1]  xa[i];
+ const double h_ip1 = xa[1]  xa[0];
const double ydiff_i = ya[i + 1]  ya[i];
 const double ydiff_ip1 = ya[(i + 2) % num_points]  ya[i + 1];
+ const double ydiff_ip1 = ya[1]  ya[0];
state>offdiag[i] = h_ip1;
state>diag[i] = 2.0 * (h_ip1 + h_i);
state>g[i] = 3.0 * (ydiff_ip1 / h_ip1  ydiff_i / h_i);
More information about the Gsldiscuss
mailing list