## File: 1D-Interpolation-Example-programs.html

package info (click to toggle)
gsl-ref-html 2.3-1
• area: non-free
• in suites: bullseye, buster, sid
• size: 6,876 kB
• ctags: 4,574
• sloc: makefile: 35
 file content (270 lines) | stat: -rw-r--r-- 9,866 bytes parent folder | download
 123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270  GNU Scientific Library – Reference Manual: 1D Interpolation Example programs

28.7 Examples of 1D Interpolation

The following program demonstrates the use of the interpolation and spline functions. It computes a cubic spline interpolation of the 10-point dataset (x_i, y_i) where x_i = i + \sin(i)/2 and y_i = i + \cos(i^2) for i = 0 \dots 9.

#include <stdlib.h> #include <stdio.h> #include <math.h> #include <gsl/gsl_errno.h> #include <gsl/gsl_spline.h>  int main (void) {   int i;   double xi, yi, x[10], y[10];    printf ("#m=0,S=2\n");    for (i = 0; i < 10; i++)     {       x[i] = i + 0.5 * sin (i);       y[i] = i + cos (i * i);       printf ("%g %g\n", x[i], y[i]);     }    printf ("#m=1,S=0\n");    {     gsl_interp_accel *acc        = gsl_interp_accel_alloc ();     gsl_spline *spline        = gsl_spline_alloc (gsl_interp_cspline, 10);      gsl_spline_init (spline, x, y, 10);      for (xi = x[0]; xi < x[9]; xi += 0.01)       {         yi = gsl_spline_eval (spline, xi, acc);         printf ("%g %g\n", xi, yi);       }     gsl_spline_free (spline);     gsl_interp_accel_free (acc);   }   return 0; }

The output is designed to be used with the GNU plotutils graph program,

$./a.out > interp.dat$ graph -T ps < interp.dat > interp.ps

The result shows a smooth interpolation of the original points. The interpolation method can be changed simply by varying the first argument of gsl_spline_alloc.

The next program demonstrates a periodic cubic spline with 4 data points. Note that the first and last points must be supplied with the same y-value for a periodic spline.

#include <stdlib.h> #include <stdio.h> #include <math.h> #include <gsl/gsl_errno.h> #include <gsl/gsl_spline.h>  int main (void) {   int N = 4;   double x[4] = {0.00, 0.10,  0.27,  0.30};   double y[4] = {0.15, 0.70, -0.10,  0.15};               /* Note: y[0] == y[3] for periodic data */    gsl_interp_accel *acc = gsl_interp_accel_alloc ();   const gsl_interp_type *t = gsl_interp_cspline_periodic;    gsl_spline *spline = gsl_spline_alloc (t, N);    int i; double xi, yi;    printf ("#m=0,S=5\n");   for (i = 0; i < N; i++)     {       printf ("%g %g\n", x[i], y[i]);     }    printf ("#m=1,S=0\n");   gsl_spline_init (spline, x, y, N);    for (i = 0; i <= 100; i++)     {       xi = (1 - i / 100.0) * x[0] + (i / 100.0) * x[N-1];       yi = gsl_spline_eval (spline, xi, acc);       printf ("%g %g\n", xi, yi);     }      gsl_spline_free (spline);   gsl_interp_accel_free (acc);   return 0; }

The output can be plotted with GNU graph.

$./a.out > interp.dat$ graph -T ps < interp.dat > interp.ps

The result shows a periodic interpolation of the original points. The slope of the fitted curve is the same at the beginning and end of the data, and the second derivative is also.

The next program illustrates the difference between the cubic spline, Akima, and Steffen interpolation types on a difficult dataset.

#include <stdio.h> #include <stdlib.h> #include <math.h>  #include <gsl/gsl_math.h> #include <gsl/gsl_spline.h>  int main(void) {   size_t i;   const size_t N = 9;    /* this dataset is taken from    * J. M. Hyman, Accurate Monotonicity preserving cubic interpolation,    * SIAM J. Sci. Stat. Comput. 4, 4, 1983. */   const double x[] = { 7.99, 8.09, 8.19, 8.7, 9.2,                        10.0, 12.0, 15.0, 20.0 };   const double y[] = { 0.0, 2.76429e-5, 4.37498e-2,                        0.169183, 0.469428, 0.943740,                        0.998636, 0.999919, 0.999994 };    gsl_interp_accel *acc = gsl_interp_accel_alloc();   gsl_spline *spline_cubic = gsl_spline_alloc(gsl_interp_cspline, N);   gsl_spline *spline_akima = gsl_spline_alloc(gsl_interp_akima, N);   gsl_spline *spline_steffen = gsl_spline_alloc(gsl_interp_steffen, N);    gsl_spline_init(spline_cubic, x, y, N);   gsl_spline_init(spline_akima, x, y, N);   gsl_spline_init(spline_steffen, x, y, N);    for (i = 0; i < N; ++i)     printf("%g %g\n", x[i], y[i]);    printf("\n\n");    for (i = 0; i <= 100; ++i)     {       double xi = (1 - i / 100.0) * x[0] + (i / 100.0) * x[N-1];       double yi_cubic = gsl_spline_eval(spline_cubic, xi, acc);       double yi_akima = gsl_spline_eval(spline_akima, xi, acc);       double yi_steffen = gsl_spline_eval(spline_steffen, xi, acc);        printf("%g %g %g %g\n", xi, yi_cubic, yi_akima, yi_steffen);     }    gsl_spline_free(spline_cubic);   gsl_spline_free(spline_akima);   gsl_spline_free(spline_steffen);   gsl_interp_accel_free(acc);    return 0; }

The cubic method exhibits a local maxima between the 6th and 7th data points and continues oscillating for the rest of the data. Akima also shows a local maxima but recovers and follows the data well after the 7th grid point. Steffen preserves monotonicity in all intervals and does not exhibit oscillations, at the expense of having a discontinuous second derivative.