| @cindex interpolation |
| @cindex spline |
| |
| 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 |
| @file{gsl_interp.h} and @file{gsl_spline.h}. |
| |
| @menu |
| * Introduction to Interpolation:: |
| * Interpolation Functions:: |
| * Interpolation Types:: |
| * Index Look-up and Acceleration:: |
| * Evaluation of Interpolating Functions:: |
| * Higher-level Interface:: |
| * Interpolation Example programs:: |
| * Interpolation References and Further Reading:: |
| @end menu |
| |
| @node Introduction to Interpolation |
| @section Introduction |
| |
| Given a set of data points @math{(x_1, y_1) \dots (x_n, y_n)} the |
| routines described in this section compute a continuous interpolating |
| function @math{y(x)} such that @math{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. |
| |
| @node Interpolation Functions |
| @section Interpolation Functions |
| |
| The interpolation function for a given dataset is stored in a |
| @code{gsl_interp} object. These are created by the following functions. |
| |
| @deftypefun {gsl_interp *} gsl_interp_alloc (const gsl_interp_type * @var{T}, size_t @var{size}) |
| This function returns a pointer to a newly allocated interpolation |
| object of type @var{T} for @var{size} data-points. |
| @end deftypefun |
| |
| @deftypefun int gsl_interp_init (gsl_interp * @var{interp}, const double @var{xa}[], const double @var{ya}[], size_t @var{size}) |
| This function initializes the interpolation object @var{interp} for the |
| data (@var{xa},@var{ya}) where @var{xa} and @var{ya} are arrays of size |
| @var{size}. The interpolation object (@code{gsl_interp}) does not save |
| the data arrays @var{xa} and @var{ya} and only stores the static state |
| computed from the data. The @var{xa} data array is always assumed to be |
| strictly ordered; the behavior for other arrangements is not defined. |
| @end deftypefun |
| |
| @deftypefun void gsl_interp_free (gsl_interp * @var{interp}) |
| This function frees the interpolation object @var{interp}. |
| @end deftypefun |
| |
| @node Interpolation Types |
| @section Interpolation Types |
| |
| The interpolation library provides five interpolation types: |
| |
| @deffn {Interpolation Type} gsl_interp_linear |
| @cindex linear interpolation |
| Linear interpolation. This interpolation method does not require any |
| additional memory. |
| @end deffn |
| |
| @deffn {Interpolation Type} gsl_interp_polynomial |
| @cindex polynomial interpolation |
| 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. |
| @end deffn |
| |
| @deffn {Interpolation Type} gsl_interp_cspline |
| @cindex cubic splines |
| 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. |
| @end deffn |
| |
| @deffn {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. |
| |
| @end deffn |
| |
| @deffn {Interpolation Type} gsl_interp_akima |
| @cindex Akima splines |
| Non-rounded Akima spline with natural boundary conditions. This method |
| uses the non-rounded corner algorithm of Wodicka. |
| @end deffn |
| |
| @deffn {Interpolation Type} gsl_interp_akima_periodic |
| Non-rounded Akima spline with periodic boundary conditions. This method |
| uses the non-rounded corner algorithm of Wodicka. |
| @end deffn |
| |
| The following related functions are available: |
| |
| @deftypefun {const char *} gsl_interp_name (const gsl_interp * @var{interp}) |
| This function returns the name of the interpolation type used by @var{interp}. |
| For example, |
| |
| @example |
| printf ("interp uses '%s' interpolation.\n", |
| gsl_interp_name (interp)); |
| @end example |
| |
| @noindent |
| would print something like, |
| |
| @example |
| interp uses 'cspline' interpolation. |
| @end example |
| @end deftypefun |
| |
| @deftypefun {unsigned int} gsl_interp_min_size (const gsl_interp * @var{interp}) |
| This function returns the minimum number of points required by the |
| interpolation type of @var{interp}. For example, Akima spline interpolation |
| requires a minimum of 5 points. |
| @end deftypefun |
| |
| @node Index Look-up and Acceleration |
| @section Index Look-up and Acceleration |
| |
| The state of searches can be stored in a @code{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. |
| |
| @deftypefun size_t gsl_interp_bsearch (const double @var{x_array}[], double @var{x}, size_t @var{index_lo}, size_t @var{index_hi}) |
| This function returns the index @math{i} of the array @var{x_array} such |
| that @code{x_array[i] <= x < x_array[i+1]}. The index is searched for |
| in the range [@var{index_lo},@var{index_hi}]. |
| @end deftypefun |
| |
| @deftypefun {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. |
| @end deftypefun |
| |
| @deftypefun size_t gsl_interp_accel_find (gsl_interp_accel * @var{a}, const double @var{x_array}[], size_t @var{size}, double @var{x}) |
| This function performs a lookup action on the data array @var{x_array} |
| of size @var{size}, using the given accelerator @var{a}. This is how |
| lookups are performed during evaluation of an interpolation. The |
| function returns an index @math{i} such that @code{x_array[i] <= x < |
| x_array[i+1]}. |
| @end deftypefun |
| |
| @deftypefun void gsl_interp_accel_free (gsl_interp_accel* @var{acc}) |
| This function frees the accelerator object @var{acc}. |
| @end deftypefun |
| |
| @node Evaluation of Interpolating Functions |
| @section Evaluation of Interpolating Functions |
| |
| @deftypefun double gsl_interp_eval (const gsl_interp * @var{interp}, const double @var{xa}[], const double @var{ya}[], double @var{x}, gsl_interp_accel * @var{acc}) |
| @deftypefunx int gsl_interp_eval_e (const gsl_interp * @var{interp}, const double @var{xa}[], const double @var{ya}[], double @var{x}, gsl_interp_accel * @var{acc}, double * @var{y}) |
| These functions return the interpolated value of @var{y} for a given |
| point @var{x}, using the interpolation object @var{interp}, data arrays |
| @var{xa} and @var{ya} and the accelerator @var{acc}. |
| @end deftypefun |
| |
| @deftypefun double gsl_interp_eval_deriv (const gsl_interp * @var{interp}, const double @var{xa}[], const double @var{ya}[], double @var{x}, gsl_interp_accel * @var{acc}) |
| @deftypefunx int gsl_interp_eval_deriv_e (const gsl_interp * @var{interp}, const double @var{xa}[], const double @var{ya}[], double @var{x}, gsl_interp_accel * @var{acc}, double * @var{d}) |
| These functions return the derivative @var{d} of an interpolated |
| function for a given point @var{x}, using the interpolation object |
| @var{interp}, data arrays @var{xa} and @var{ya} and the accelerator |
| @var{acc}. |
| @end deftypefun |
| |
| @deftypefun double gsl_interp_eval_deriv2 (const gsl_interp * @var{interp}, const double @var{xa}[], const double @var{ya}[], double @var{x}, gsl_interp_accel * @var{acc}) |
| @deftypefunx int gsl_interp_eval_deriv2_e (const gsl_interp * @var{interp}, const double @var{xa}[], const double @var{ya}[], double @var{x}, gsl_interp_accel * @var{acc}, double * @var{d2}) |
| These functions return the second derivative @var{d2} of an interpolated |
| function for a given point @var{x}, using the interpolation object |
| @var{interp}, data arrays @var{xa} and @var{ya} and the accelerator |
| @var{acc}. |
| @end deftypefun |
| |
| @deftypefun double gsl_interp_eval_integ (const gsl_interp * @var{interp}, const double @var{xa}[], const double @var{ya}[], double @var{a}, double @var{b}, gsl_interp_accel * @var{acc}) |
| @deftypefunx int gsl_interp_eval_integ_e (const gsl_interp * @var{interp}, const double @var{xa}[], const double @var{ya}[], double @var{a}, double @var{b}, gsl_interp_accel * @var{acc}, double * @var{result}) |
| These functions return the numerical integral @var{result} of an |
| interpolated function over the range [@var{a}, @var{b}], using the |
| interpolation object @var{interp}, data arrays @var{xa} and @var{ya} and |
| the accelerator @var{acc}. |
| @end deftypefun |
| |
| @node Higher-level Interface |
| @section Higher-level Interface |
| |
| The functions described in the previous sections required the user to |
| supply pointers to the @math{x} and @math{y} arrays on each call. The |
| following functions are equivalent to the corresponding |
| @code{gsl_interp} functions but maintain a copy of this data in the |
| @code{gsl_spline} object. This removes the need to pass both @var{xa} |
| and @var{ya} as arguments on each evaluation. These functions are |
| defined in the header file @file{gsl_spline.h}. |
| |
| @deftypefun {gsl_spline *} gsl_spline_alloc (const gsl_interp_type * @var{T}, size_t @var{size}) |
| @end deftypefun |
| |
| @deftypefun int gsl_spline_init (gsl_spline * @var{spline}, const double @var{xa}[], const double @var{ya}[], size_t @var{size}) |
| @end deftypefun |
| |
| @deftypefun void gsl_spline_free (gsl_spline * @var{spline}) |
| @end deftypefun |
| |
| @deftypefun {const char *} gsl_spline_name (const gsl_spline * @var{spline}) |
| @end deftypefun |
| |
| @deftypefun {unsigned int} gsl_spline_min_size (const gsl_spline * @var{spline}) |
| @end deftypefun |
| |
| @deftypefun double gsl_spline_eval (const gsl_spline * @var{spline}, double @var{x}, gsl_interp_accel * @var{acc}) |
| @deftypefunx int gsl_spline_eval_e (const gsl_spline * @var{spline}, double @var{x}, gsl_interp_accel * @var{acc}, double * @var{y}) |
| @end deftypefun |
| |
| @deftypefun double gsl_spline_eval_deriv (const gsl_spline * @var{spline}, double @var{x}, gsl_interp_accel * @var{acc}) |
| @deftypefunx int gsl_spline_eval_deriv_e (const gsl_spline * @var{spline}, double @var{x}, gsl_interp_accel * @var{acc}, double * @var{d}) |
| @end deftypefun |
| |
| @deftypefun double gsl_spline_eval_deriv2 (const gsl_spline * @var{spline}, double @var{x}, gsl_interp_accel * @var{acc}) |
| @deftypefunx int gsl_spline_eval_deriv2_e (const gsl_spline * @var{spline}, double @var{x}, gsl_interp_accel * @var{acc}, double * @var{d2}) |
| @end deftypefun |
| |
| @deftypefun double gsl_spline_eval_integ (const gsl_spline * @var{spline}, double @var{a}, double @var{b}, gsl_interp_accel * @var{acc}) |
| @deftypefunx int gsl_spline_eval_integ_e (const gsl_spline * @var{spline}, double @var{a}, double @var{b}, gsl_interp_accel * @var{acc}, double * @var{result}) |
| @end deftypefun |
| |
| @node Interpolation Example programs |
| @section Examples |
| |
| The following program demonstrates the use of the interpolation and |
| spline functions. It computes a cubic spline interpolation of the |
| 10-point dataset @math{(x_i, y_i)} where @math{x_i = i + \sin(i)/2} and |
| @math{y_i = i + \cos(i^2)} for @math{i = 0 \dots 9}. |
| |
| @example |
| @verbatiminclude examples/interp.c |
| @end example |
| |
| @noindent |
| The output is designed to be used with the @sc{gnu} plotutils |
| @code{graph} program, |
| |
| @example |
| $ ./a.out > interp.dat |
| $ graph -T ps < interp.dat > interp.ps |
| @end example |
| |
| @iftex |
| @sp 1 |
| @center @image{interp2,3.4in} |
| @end iftex |
| |
| @noindent |
| The result shows a smooth interpolation of the original points. The |
| interpolation method can changed simply by varying the first argument of |
| @code{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. |
| |
| @example |
| @verbatiminclude examples/interpp.c |
| @end example |
| |
| @noindent |
| |
| The output can be plotted with @sc{gnu} @code{graph}. |
| |
| @example |
| $ ./a.out > interp.dat |
| $ graph -T ps < interp.dat > interp.ps |
| @end example |
| |
| @iftex |
| @sp 1 |
| @center @image{interpp2,3.4in} |
| @end iftex |
| |
| @noindent |
| 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. |
| |
| @node Interpolation References and Further Reading |
| @section References and Further Reading |
| |
| Descriptions of the interpolation algorithms and further references can |
| be found in the following books: |
| |
| @itemize @asis |
| @item C.W. Ueberhuber, |
| @cite{Numerical Computation (Volume 1), Chapter 9 ``Interpolation''}, |
| Springer (1997), ISBN 3-540-62058-3. |
| |
| @item D.M. Young, R.T. Gregory |
| @cite{A Survey of Numerical Mathematics (Volume 1), Chapter 6.8}, |
| Dover (1988), ISBN 0-486-65691-8. |
| @end itemize |
| |
| @noindent |