htswanalysis/MACS/lib/gsl/gsl-1.11/doc/interp.texi

   1 @cindex interpolation
   2 @cindex spline
   3
   4 This chapter describes functions for performing interpolation.  The
   5 library provides a variety of interpolation methods, including Cubic
   6 splines and Akima splines.  The interpolation types are interchangeable,
   7 allowing different methods to be used without recompiling.
   8 Interpolations can be defined for both normal and periodic boundary
   9 conditions.  Additional functions are available for computing
  10 derivatives and integrals of interpolating functions.
  11
  12 The functions described in this section are declared in the header files
  13 @file{gsl_interp.h} and @file{gsl_spline.h}.
  14
  15 @menu
  16 * Introduction to Interpolation::
  17 * Interpolation Functions::
  18 * Interpolation Types::
  19 * Index Look-up and Acceleration::
  20 * Evaluation of Interpolating Functions::
  21 * Higher-level Interface::
  22 * Interpolation Example programs::
  23 * Interpolation References and Further Reading::
  24 @end menu
  25
  26 @node Introduction to Interpolation
  27 @section Introduction
  28
  29 Given a set of data points @math{(x_1, y_1) \dots (x_n, y_n)} the
  30 routines described in this section compute a continuous interpolating
  31 function @math{y(x)} such that @math{y(x_i) = y_i}.  The interpolation
  32 is piecewise smooth, and its behavior at the end-points is determined by
  33 the type of interpolation used.
  34
  35 @node Interpolation Functions
  36 @section Interpolation Functions
  37
  38 The interpolation function for a given dataset is stored in a
  39 @code{gsl_interp} object.  These are created by the following functions.
  40
  41 @deftypefun {gsl_interp *} gsl_interp_alloc (const gsl_interp_type * @var{T}, size_t @var{size})
  42 This function returns a pointer to a newly allocated interpolation
  43 object of type @var{T} for @var{size} data-points.
  44 @end deftypefun
  45
  46 @deftypefun int gsl_interp_init (gsl_interp * @var{interp}, const double @var{xa}[], const double @var{ya}[], size_t @var{size})
  47 This function initializes the interpolation object @var{interp} for the
  48 data (@var{xa},@var{ya}) where @var{xa} and @var{ya} are arrays of size
  49 @var{size}.  The interpolation object (@code{gsl_interp}) does not save
  50 the data arrays @var{xa} and @var{ya} and only stores the static state
  51 computed from the data.  The @var{xa} data array is always assumed to be
  52 strictly ordered, with increasing @math{x} values;
  53 the behavior for other arrangements is not defined.
  54 @end deftypefun
  55
  56 @deftypefun void gsl_interp_free (gsl_interp * @var{interp})
  57 This function frees the interpolation object @var{interp}.
  58 @end deftypefun
  59
  60 @node Interpolation Types
  61 @section Interpolation Types
  62
  63 The interpolation library provides five interpolation types:
  64
  65 @deffn {Interpolation Type} gsl_interp_linear
  66 @cindex linear interpolation
  67 Linear interpolation.  This interpolation method does not require any
  68 additional memory.
  69 @end deffn
  70
  71 @deffn {Interpolation Type} gsl_interp_polynomial
  72 @cindex polynomial interpolation
  73 Polynomial interpolation.  This method should only be used for
  74 interpolating small numbers of points because polynomial interpolation
  75 introduces large oscillations, even for well-behaved datasets.  The
  76 number of terms in the interpolating polynomial is equal to the number
  77 of points.
  78 @end deffn
  79
  80 @deffn {Interpolation Type} gsl_interp_cspline
  81 @cindex cubic splines
  82 Cubic spline with natural boundary conditions.  The resulting curve is
  83 piecewise cubic on each interval, with matching first and second
  84 derivatives at the supplied data-points.  The second derivative is
  85 chosen to be zero at the first point and last point.
  86 @end deffn
  87
  88 @deffn {Interpolation Type} gsl_interp_cspline_periodic
  89 Cubic spline with periodic boundary conditions.  The resulting curve
  90 is piecewise cubic on each interval, with matching first and second
  91 derivatives at the supplied data-points.  The derivatives at the first
  92 and last points are also matched.  Note that the last point in the
  93 data must have the same y-value as the first point, otherwise the
  94 resulting periodic interpolation will have a discontinuity at the
  95 boundary.
  96
  97 @end deffn
  98
  99 @deffn {Interpolation Type} gsl_interp_akima
 100 @cindex Akima splines
 101 Non-rounded Akima spline with natural boundary conditions.  This method
 102 uses the non-rounded corner algorithm of Wodicka.
 103 @end deffn
 104
 105 @deffn {Interpolation Type} gsl_interp_akima_periodic
 106 Non-rounded Akima spline with periodic boundary conditions.  This method
 107 uses the non-rounded corner algorithm of Wodicka.
 108 @end deffn
 109
 110 The following related functions are available:
 111
 112 @deftypefun {const char *} gsl_interp_name (const gsl_interp * @var{interp})
 113 This function returns the name of the interpolation type used by @var{interp}.
 114 For example,
 115
 116 @example
 117 printf ("interp uses '%s' interpolation.\n",
 118         gsl_interp_name (interp));
 119 @end example
 120
 121 @noindent
 122 would print something like,
 123
 124 @example
 125 interp uses 'cspline' interpolation.
 126 @end example
 127 @end deftypefun
 128
 129 @deftypefun {unsigned int} gsl_interp_min_size (const gsl_interp * @var{interp})
 130 This function returns the minimum number of points required by the
 131 interpolation type of @var{interp}.  For example, Akima spline interpolation
 132 requires a minimum of 5 points.
 133 @end deftypefun
 134
 135 @node Index Look-up and Acceleration
 136 @section Index Look-up and Acceleration
 137
 138 The state of searches can be stored in a @code{gsl_interp_accel} object,
 139 which is a kind of iterator for interpolation lookups.  It caches the
 140 previous value of an index lookup.  When the subsequent interpolation
 141 point falls in the same interval its index value can be returned
 142 immediately.
 143
 144 @deftypefun size_t gsl_interp_bsearch (const double @var{x_array}[], double @var{x}, size_t @var{index_lo}, size_t @var{index_hi})
 145 This function returns the index @math{i} of the array @var{x_array} such
 146 that @code{x_array[i] <= x < x_array[i+1]}.  The index is searched for
 147 in the range [@var{index_lo},@var{index_hi}].  @inlinefn{}
 148 @end deftypefun
 149
 150 @deftypefun {gsl_interp_accel *} gsl_interp_accel_alloc (void)
 151 This function returns a pointer to an accelerator object, which is a
 152 kind of iterator for interpolation lookups.  It tracks the state of
 153 lookups, thus allowing for application of various acceleration
 154 strategies.
 155 @end deftypefun
 156
 157 @deftypefun size_t gsl_interp_accel_find (gsl_interp_accel * @var{a}, const double @var{x_array}[], size_t @var{size}, double @var{x})
 158 This function performs a lookup action on the data array @var{x_array}
 159 of size @var{size}, using the given accelerator @var{a}.  This is how
 160 lookups are performed during evaluation of an interpolation.  The
 161 function returns an index @math{i} such that @code{x_array[i] <= x <
 162 x_array[i+1]}.  @inlinefn{}
 163 @end deftypefun
 164
 165 @deftypefun void gsl_interp_accel_free (gsl_interp_accel* @var{acc})
 166 This function frees the accelerator object @var{acc}.
 167 @end deftypefun
 168
 169 @node Evaluation of Interpolating Functions
 170 @section Evaluation of Interpolating Functions
 171
 172 @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})
 173 @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})
 174 These functions return the interpolated value of @var{y} for a given
 175 point @var{x}, using the interpolation object @var{interp}, data arrays
 176 @var{xa} and @var{ya} and the accelerator @var{acc}.
 177 @end deftypefun
 178
 179 @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})
 180 @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})
 181 These functions return the derivative @var{d} of an interpolated
 182 function for a given point @var{x}, using the interpolation object
 183 @var{interp}, data arrays @var{xa} and @var{ya} and the accelerator
 184 @var{acc}.
 185 @end deftypefun
 186
 187 @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})
 188 @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})
 189 These functions return the second derivative @var{d2} of an interpolated
 190 function for a given point @var{x}, using the interpolation object
 191 @var{interp}, data arrays @var{xa} and @var{ya} and the accelerator
 192 @var{acc}.
 193 @end deftypefun
 194
 195 @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})
 196 @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})
 197 These functions return the numerical integral @var{result} of an
 198 interpolated function over the range [@var{a}, @var{b}], using the
 199 interpolation object @var{interp}, data arrays @var{xa} and @var{ya} and
 200 the accelerator @var{acc}.
 201 @end deftypefun
 202
 203 @node Higher-level Interface
 204 @section Higher-level Interface
 205
 206 The functions described in the previous sections required the user to
 207 supply pointers to the @math{x} and @math{y} arrays on each call.  The
 208 following functions are equivalent to the corresponding
 209 @code{gsl_interp} functions but maintain a copy of this data in the
 210 @code{gsl_spline} object.  This removes the need to pass both @var{xa}
 211 and @var{ya} as arguments on each evaluation. These functions are
 212 defined in the header file @file{gsl_spline.h}.
 213
 214 @deftypefun {gsl_spline *} gsl_spline_alloc (const gsl_interp_type * @var{T}, size_t @var{size})
 215 @end deftypefun
 216
 217 @deftypefun int gsl_spline_init (gsl_spline * @var{spline}, const double @var{xa}[], const double @var{ya}[], size_t @var{size})
 218 @end deftypefun
 219
 220 @deftypefun void gsl_spline_free (gsl_spline * @var{spline})
 221 @end deftypefun
 222
 223 @deftypefun {const char *} gsl_spline_name (const gsl_spline * @var{spline})
 224 @end deftypefun
 225
 226 @deftypefun {unsigned int} gsl_spline_min_size (const gsl_spline * @var{spline})
 227 @end deftypefun
 228
 229 @deftypefun double gsl_spline_eval (const gsl_spline * @var{spline}, double @var{x}, gsl_interp_accel * @var{acc})
 230 @deftypefunx int gsl_spline_eval_e (const gsl_spline * @var{spline}, double @var{x}, gsl_interp_accel * @var{acc}, double * @var{y})
 231 @end deftypefun
 232
 233 @deftypefun double gsl_spline_eval_deriv (const gsl_spline * @var{spline}, double @var{x}, gsl_interp_accel * @var{acc})
 234 @deftypefunx int gsl_spline_eval_deriv_e (const gsl_spline * @var{spline}, double @var{x}, gsl_interp_accel * @var{acc}, double * @var{d})
 235 @end deftypefun
 236
 237 @deftypefun double gsl_spline_eval_deriv2 (const gsl_spline * @var{spline}, double @var{x}, gsl_interp_accel * @var{acc})
 238 @deftypefunx int gsl_spline_eval_deriv2_e (const gsl_spline * @var{spline}, double @var{x}, gsl_interp_accel * @var{acc}, double * @var{d2})
 239 @end deftypefun
 240
 241 @deftypefun double gsl_spline_eval_integ (const gsl_spline * @var{spline}, double @var{a}, double @var{b}, gsl_interp_accel * @var{acc})
 242 @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})
 243 @end deftypefun
 244
 245 @node Interpolation Example programs
 246 @section Examples
 247
 248 The following program demonstrates the use of the interpolation and
 249 spline functions.  It computes a cubic spline interpolation of the
 250 10-point dataset @math{(x_i, y_i)} where @math{x_i = i + \sin(i)/2} and
 251 @math{y_i = i + \cos(i^2)} for @math{i = 0 \dots 9}.
 252
 253 @example
 254 @verbatiminclude examples/interp.c
 255 @end example
 256
 257 @noindent
 258 The output is designed to be used with the @sc{gnu} plotutils
 259 @code{graph} program,
 260
 261 @example
 262 $ ./a.out > interp.dat
 263 $ graph -T ps < interp.dat > interp.ps
 264 @end example
 265
 266 @iftex
 267 @sp 1
 268 @center @image{interp2,3.4in}
 269 @end iftex
 270
 271 @noindent
 272 The result shows a smooth interpolation of the original points.  The
 273 interpolation method can be changed simply by varying the first argument of
 274 @code{gsl_spline_alloc}.
 275
 276 The next program demonstrates a periodic cubic spline with 4 data
 277 points.  Note that the first and last points must be supplied with
 278 the same y-value for a periodic spline.
 279
 280 @example
 281 @verbatiminclude examples/interpp.c
 282 @end example
 283
 284 @noindent
 285
 286 The output can be plotted with @sc{gnu} @code{graph}.
 287
 288 @example
 289 $ ./a.out > interp.dat
 290 $ graph -T ps < interp.dat > interp.ps
 291 @end example
 292
 293 @iftex
 294 @sp 1
 295 @center @image{interpp2,3.4in}
 296 @end iftex
 297
 298 @noindent
 299 The result shows a periodic interpolation of the original points. The
 300 slope of the fitted curve is the same at the beginning and end of the
 301 data, and the second derivative is also.
 302
 303 @node Interpolation References and Further Reading
 304 @section References and Further Reading
 305
 306 Descriptions of the interpolation algorithms and further references can
 307 be found in the following books:
 308
 309 @itemize @asis
 310 @item C.W. Ueberhuber,
 311 @cite{Numerical Computation (Volume 1), Chapter 9 ``Interpolation''},
 312 Springer (1997), ISBN 3-540-62058-3.
 313
 314 @item D.M. Young, R.T. Gregory
 315 @cite{A Survey of Numerical Mathematics (Volume 1), Chapter 6.8},
 316 Dover (1988), ISBN 0-486-65691-8.
 317 @end itemize
 318
 319 @noindent