This chapter describes functions for performing interpolation. The
library provides a variety of interpolation methods, including Cubic
splines and Akima splines. The interpolation types are interchangeable,
allowing different methods to be used without recompiling.
Interpolations can be defined for both normal and periodic boundary
conditions. Additional functions are available for computing
derivatives and integrals of interpolating functions.
The functions described in this section are declared in the header files
`gsl_interp.h' and `gsl_spline.h'.
Given a set of data points (x_1, y_1) \dots (x_n, y_n) the
routines described in this section compute a continuous interpolating
function y(x) such that y_i = y(x_i). The interpolation
is piecewise smooth, and its behavior at the points is determined by the
type of interpolation used.
This function initializes the interpolation object interp for the
data (xa,ya) where xa and ya are arrays of size
size. The interpolation object (gsl_interp) does not save
the data arrays xa and ya and only stores the static state
computed from the data. The xa data array is always assumed to be
strictly ordered; the behavior for other arrangements is not defined.
The interpolation library provides five interpolation types:
Interpolation Type:gsl_interp_linear
Linear interpolation. This interpolation method does not require any
additional memory.
Interpolation Type:gsl_interp_polynomial
Polynomial interpolation. This method should only be used for
interpolating small numbers of points because polynomial interpolation
introduces large oscillations, even for well-behaved datasets. The
number of terms in the interpolating polynomial is equal to the number
of points.
Function: unsigned int gsl_interp_min_size(const gsl_interp * interp)
This function returns the minimum number of points required by the
interpolation type of interp. For example, Akima spline interpolation
requires a minimum of 5 points.
The state of searches can be stored in a gsl_interp_accel object,
which is a kind of iterator for interpolation lookups. It caches the
previous value of an index lookup. When the subsequent interpolation
point falls in the same interval its index value can be returned
immediately.
This function returns the index i of the array x_array such
that x_array[i] <= x < x_array[i+1]. The index is searched for
in the range [index_lo,index_hi].
This function returns a pointer to an accelerator object, which is a
kind of iterator for interpolation lookups. It tracks the state of
lookups, thus allowing for application of various acceleration
strategies.
This function performs a lookup action on the data array x_array
of size size, using the given accelerator a. This is how
lookups are performed during evaluation of an interpolation. The
function returns an index i such that xarray[i] <= x <
xarray[i+1].
Function: void gsl_interp_accel_free(gsl_interp_accel* a)
Function: int gsl_interp_eval_e(const gsl_interp * interp, const double xa[], const double ya[], double x, gsl_interp_accel * a, double * y)
These functions return the interpolated value of y for a given
point x, using the interpolation object interp, data arrays
xa and ya and the accelerator a.
Function: int gsl_interp_eval_deriv_e(const gsl_interp * interp, const double xa[], const double ya[], double x, gsl_interp_accel * a, double * d)
These functions return the derivative d of an interpolated
function for a given point x, using the interpolation object
interp, data arrays xa and ya and the accelerator
a.
Function: int gsl_interp_eval_deriv2_e(const gsl_interp * interp, const double xa[], const double ya[], double x, gsl_interp_accel * a, double * d2)
These functions return the second derivative d2 of an interpolated
function for a given point x, using the interpolation object
interp, data arrays xa and ya and the accelerator
a.
Function: double gsl_interp_eval_integ(const gsl_interp * interp, const double xa[], const double ya[], double a, double bdouble x, gsl_interp_accel * a)
Function: int gsl_interp_eval_integ_e(const gsl_interp * interp, const double xa[], const double ya[], , double a, double b, gsl_interp_accel * a, double * result)
These functions return the numerical integral result of an
interpolated function over the range [a, b], using the
interpolation object interp, data arrays xa and ya and
the accelerator a.
The functions described in the previous sections required the user to
supply pointers to the x and y arrays on each call. The
following functions are equivalent to the corresponding
gsl_interp functions but maintain a copy of this data in the
gsl_spline object. This removes the need to pass both xa
and ya as arguments on each evaluation. These functions are
defined in the header file `gsl_spline.h'.
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 <config.h>
#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,
The result shows a smooth interpolation of the original points. The
interpolation method can changed simply by varying the first argument of
gsl_spline_alloc.