Go to the first, previous, next, last section, table of contents.

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.

These interpolation methods produce curves that pass through each datapoint. To interpolate noisy data with a smoothing curve see section Basis Splines.

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(x_i) = y_i. The interpolation is piecewise smooth, and its behavior at the end-points is determined by the type of interpolation used.

The interpolation function for a given dataset is stored in a
`gsl_interp`

object. These are created by the following functions.

__Function:__gsl_interp ***gsl_interp_alloc***(const gsl_interp_type **`T`, size_t`size`)-
This function returns a pointer to a newly allocated interpolation
object of type
`T`for`size`data-points.

__Function:__int**gsl_interp_init***(gsl_interp **`interp`, const double`xa`[], const double`ya`[], size_t`size`)-
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, with increasing x values; the behavior for other arrangements is not defined.

__Function:__void**gsl_interp_free***(gsl_interp **`interp`)-
This function frees the interpolation object
`interp`.

The interpolation library provides six 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.

__Interpolation Type:__**gsl_interp_cspline**- Cubic spline with natural boundary conditions. The resulting curve is piecewise cubic on each interval, with matching first and second derivatives at the supplied data-points. The second derivative is chosen to be zero at the first point and last point.

__Interpolation Type:__**gsl_interp_cspline_periodic**- Cubic spline with periodic boundary conditions. The resulting curve is piecewise cubic on each interval, with matching first and second derivatives at the supplied data-points. The derivatives at the first and last points are also matched. Note that the last point in the data must have the same y-value as the first point, otherwise the resulting periodic interpolation will have a discontinuity at the boundary.

__Interpolation Type:__**gsl_interp_akima**- Non-rounded Akima spline with natural boundary conditions. This method uses the non-rounded corner algorithm of Wodicka.

__Interpolation Type:__**gsl_interp_akima_periodic**- Non-rounded Akima spline with periodic boundary conditions. This method uses the non-rounded corner algorithm of Wodicka.

The following related functions are available:

__Function:__const char ***gsl_interp_name***(const gsl_interp **`interp`)-
This function returns the name of the interpolation type used by
`interp`. For example,printf ("interp uses '%s' interpolation.\n", gsl_interp_name (interp));

would print something like,

interp uses 'cspline' interpolation.

__Function:__unsigned int**gsl_interp_min_size***(const gsl_interp **`interp`)__Function:__unsigned int**gsl_interp_type_min_size***(const gsl_interp_type **`T`)-
These functions return the minimum number of points required by the
interpolation object
`interp`or interpolation type`T`. 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.

__Function:__size_t**gsl_interp_bsearch***(const double*`x_array`[], double`x`, size_t`index_lo`, size_t`index_hi`)-
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`]. @inlinefn{}

__Function:__gsl_interp_accel ***gsl_interp_accel_alloc***(void)*- 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.

__Function:__size_t**gsl_interp_accel_find***(gsl_interp_accel **`a`, const double`x_array`[], size_t`size`, double`x`)-
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`x_array[i] <= x < x_array[i+1]`

. @inlinefn{}

__Function:__int**gsl_interp_accel_reset***(gsl_interp_accel **`acc`);-
This function reinitializes the accelerator object
`acc`. It should be used when the cached information is no longer applicable--for example, when switching to a new dataset.

__Function:__void**gsl_interp_accel_free***(gsl_interp_accel**`acc`)-
This function frees the accelerator object
`acc`.

__Function:__double**gsl_interp_eval***(const gsl_interp **`interp`, const double`xa`[], const double`ya`[], double`x`, gsl_interp_accel *`acc`)__Function:__int**gsl_interp_eval_e***(const gsl_interp **`interp`, const double`xa`[], const double`ya`[], double`x`, gsl_interp_accel *`acc`, 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`acc`. When`x`is outside the range of`xa`, the error code`GSL_EDOM`

is returned with a value of`GSL_NAN`

for`y`.

__Function:__double**gsl_interp_eval_deriv***(const gsl_interp **`interp`, const double`xa`[], const double`ya`[], double`x`, gsl_interp_accel *`acc`)__Function:__int**gsl_interp_eval_deriv_e***(const gsl_interp **`interp`, const double`xa`[], const double`ya`[], double`x`, gsl_interp_accel *`acc`, 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`acc`.

__Function:__double**gsl_interp_eval_deriv2***(const gsl_interp **`interp`, const double`xa`[], const double`ya`[], double`x`, gsl_interp_accel *`acc`)__Function:__int**gsl_interp_eval_deriv2_e***(const gsl_interp **`interp`, const double`xa`[], const double`ya`[], double`x`, gsl_interp_accel *`acc`, 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`acc`.

__Function:__double**gsl_interp_eval_integ***(const gsl_interp **`interp`, const double`xa`[], const double`ya`[], double`a`, double`b`, gsl_interp_accel *`acc`)__Function:__int**gsl_interp_eval_integ_e***(const gsl_interp **`interp`, const double`xa`[], const double`ya`[], double`a`, double`b`, gsl_interp_accel *`acc`, 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`acc`.

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'`.

__Function:__gsl_spline ***gsl_spline_alloc***(const gsl_interp_type **`T`, size_t`size`)

__Function:__int**gsl_spline_init***(gsl_spline **`spline`, const double`xa`[], const double`ya`[], size_t`size`)

__Function:__void**gsl_spline_free***(gsl_spline **`spline`)

__Function:__const char ***gsl_spline_name***(const gsl_spline **`spline`)

__Function:__unsigned int**gsl_spline_min_size***(const gsl_spline **`spline`)

__Function:__double**gsl_spline_eval***(const gsl_spline **`spline`, double`x`, gsl_interp_accel *`acc`)__Function:__int**gsl_spline_eval_e***(const gsl_spline **`spline`, double`x`, gsl_interp_accel *`acc`, double *`y`)

__Function:__double**gsl_spline_eval_deriv***(const gsl_spline **`spline`, double`x`, gsl_interp_accel *`acc`)__Function:__int**gsl_spline_eval_deriv_e***(const gsl_spline **`spline`, double`x`, gsl_interp_accel *`acc`, double *`d`)

__Function:__double**gsl_spline_eval_deriv2***(const gsl_spline **`spline`, double`x`, gsl_interp_accel *`acc`)__Function:__int**gsl_spline_eval_deriv2_e***(const gsl_spline **`spline`, double`x`, gsl_interp_accel *`acc`, double *`d2`)

__Function:__double**gsl_spline_eval_integ***(const gsl_spline **`spline`, double`a`, double`b`, gsl_interp_accel *`acc`)__Function:__int**gsl_spline_eval_integ_e***(const gsl_spline **`spline`, double`a`, double`b`, gsl_interp_accel *`acc`, double *`result`)

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.

Descriptions of the interpolation algorithms and further references can be found in the following books:

- C.W. Ueberhuber, Numerical Computation (Volume 1), Chapter 9 "Interpolation", Springer (1997), ISBN 3-540-62058-3.
- D.M. Young, R.T. Gregory A Survey of Numerical Mathematics (Volume 1), Chapter 6.8, Dover (1988), ISBN 0-486-65691-8.

Go to the first, previous, next, last section, table of contents.