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

Interpolation

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

Introduction

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 end-points is determined by the type of interpolation used.

Interpolation Functions

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.

Interpolation Types

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: 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): 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.

Index Look-up and Acceleration

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 x_array[i] <= x < x_array[i+1].

Function: void gsl_interp_accel_free (gsl_interp_accel* acc): This function frees the accelerator object acc.

Evaluation of Interpolating Functions

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.

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.

Higher-level Interface

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: 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)

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 changed simply by varying the first argument of gsl_spline_alloc.

References and Further Reading

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.