26.7 Examples
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: first = last
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.
