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.
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.
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; 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 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.
- Interpolation Type: gsl_interp_cspline
-
Cubic spline with natural boundary conditions
- Interpolation Type: gsl_interp_cspline_periodic
-
Cubic spline with periodic boundary conditions
- Interpolation Type: gsl_interp_akima
-
Akima spline with natural boundary conditions
- Interpolation Type: gsl_interp_akima_periodic
-
Akima spline with periodic boundary conditions
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)
-
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.
- 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].
- 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
xarray[i] <= x <
xarray[i+1]
.
- Function: void gsl_interp_accel_free (gsl_interp_accel* a)
-
This function frees the accelerator object a.
- Function: double gsl_interp_eval (const gsl_interp * interp, const double xa[], const double ya[], double x, 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: double gsl_interp_eval_deriv (const gsl_interp * interp, const double xa[], const double ya[], double x, gsl_interp_accel * 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: double gsl_interp_eval_deriv2 (const gsl_interp * interp, const double xa[], const double ya[], double x, gsl_interp_accel * 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'.
- Function: gsl_spline * gsl_spline_alloc (const gsl_interp_type * T, size_t n)
-
- 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: double gsl_spline_eval (const gsl_spline * spline, double x, gsl_interp_accel * a)
-
- Function: int gsl_spline_eval_e (const gsl_spline * spline, double x, gsl_interp_accel * a, double * y)
-
- Function: double gsl_spline_eval_deriv (const gsl_spline * spline, double x, gsl_interp_accel * a)
-
- Function: int gsl_spline_eval_deriv_e (const gsl_spline * spline, double x, gsl_interp_accel * a, double * d)
-
- Function: double gsl_spline_eval_deriv2 (const gsl_spline * spline, double x, gsl_interp_accel * a)
-
- Function: int gsl_spline_eval_deriv2_e (const gsl_spline * spline, double x, gsl_interp_accel * a, 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 <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,
$ ./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 changed simply by varying the first argument of
gsl_spline_alloc
.
Descriptions of the interpolation algorithms and further references can
be found in the following book,
- C.W. Ueberhuber,
Numerical Computation (Volume 1), Chapter 9 "Interpolation",
Springer (1997), ISBN 3-540-62058-3.
Go to the first, previous, next, last section, table of contents.