~ubuntu-branches/ubuntu/dapper/gsl-ref-html/dapper : contents of gsl-ref

~ubuntu-branches/ubuntu/dapper/gsl-ref-html/dapper : (revision 4)
<HTML>
<HEAD>
<!-- This HTML file has been created by texi2html 1.54+ (gsl)
     from /home/bjg/gsl.redhat/doc/gsl-ref.texi on 14 September 2005 -->

<TITLE>GNU Scientific Library -- Reference Manual - Interpolation</TITLE>
<link href="gsl-ref_27.html" rel=Next>
<link href="gsl-ref_25.html" rel=Previous>
<link href="gsl-ref_toc.html" rel=ToC>

</HEAD>
<BODY>
<p>Go to the <A HREF="gsl-ref_1.html">first</A>, <A HREF="gsl-ref_25.html">previous</A>, <A HREF="gsl-ref_27.html">next</A>, <A HREF="gsl-ref_50.html">last</A> section, <A HREF="gsl-ref_toc.html">table of contents</A>.
<P><HR><P>


<H1><A NAME="SEC398" HREF="gsl-ref_toc.html#TOC398">Interpolation</A></H1>
<P>
<A NAME="IDX1945"></A>
<A NAME="IDX1946"></A>

</P>
<P>
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.

</P>
<P>
The functions described in this section are declared in the header files
<TT>`gsl_interp.h'</TT> and <TT>`gsl_spline.h'</TT>.

</P>



<H2><A NAME="SEC399" HREF="gsl-ref_toc.html#TOC399">Introduction</A></H2>

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

</P>


<H2><A NAME="SEC400" HREF="gsl-ref_toc.html#TOC400">Interpolation Functions</A></H2>

<P>
The interpolation function for a given dataset is stored in a
<CODE>gsl_interp</CODE> object.  These are created by the following functions.

</P>
<P>
<DL>
<DT><U>Function:</U> gsl_interp * <B>gsl_interp_alloc</B> <I>(const gsl_interp_type * <VAR>T</VAR>, size_t <VAR>size</VAR>)</I>
<DD><A NAME="IDX1947"></A>
This function returns a pointer to a newly allocated interpolation
object of type <VAR>T</VAR> for <VAR>size</VAR> data-points.
</DL>

</P>
<P>
<DL>
<DT><U>Function:</U> int <B>gsl_interp_init</B> <I>(gsl_interp * <VAR>interp</VAR>, const double <VAR>xa</VAR>[], const double <VAR>ya</VAR>[], size_t <VAR>size</VAR>)</I>
<DD><A NAME="IDX1948"></A>
This function initializes the interpolation object <VAR>interp</VAR> for the
data (<VAR>xa</VAR>,<VAR>ya</VAR>) where <VAR>xa</VAR> and <VAR>ya</VAR> are arrays of size
<VAR>size</VAR>.  The interpolation object (<CODE>gsl_interp</CODE>) does not save
the data arrays <VAR>xa</VAR> and <VAR>ya</VAR> and only stores the static state
computed from the data.  The <VAR>xa</VAR> data array is always assumed to be
strictly ordered; the behavior for other arrangements is not defined.
</DL>

</P>
<P>
<DL>
<DT><U>Function:</U> void <B>gsl_interp_free</B> <I>(gsl_interp * <VAR>interp</VAR>)</I>
<DD><A NAME="IDX1949"></A>
This function frees the interpolation object <VAR>interp</VAR>.
</DL>

</P>


<H2><A NAME="SEC401" HREF="gsl-ref_toc.html#TOC401">Interpolation Types</A></H2>

<P>
The interpolation library provides five interpolation types:

</P>
<P>
<DL>
<DT><U>Interpolation Type:</U> <B>gsl_interp_linear</B>
<DD><A NAME="IDX1950"></A>
<A NAME="IDX1951"></A>
Linear interpolation.  This interpolation method does not require any
additional memory.
</DL>

</P>
<P>
<DL>
<DT><U>Interpolation Type:</U> <B>gsl_interp_polynomial</B>
<DD><A NAME="IDX1952"></A>
<A NAME="IDX1953"></A>
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.
</DL>

</P>
<P>
<DL>
<DT><U>Interpolation Type:</U> <B>gsl_interp_cspline</B>
<DD><A NAME="IDX1954"></A>
<A NAME="IDX1955"></A>
Cubic spline with natural boundary conditions.
</DL>

</P>
<P>
<DL>
<DT><U>Interpolation Type:</U> <B>gsl_interp_cspline_periodic</B>
<DD><A NAME="IDX1956"></A>
Cubic spline with periodic boundary conditions
</DL>

</P>
<P>
<DL>
<DT><U>Interpolation Type:</U> <B>gsl_interp_akima</B>
<DD><A NAME="IDX1957"></A>
<A NAME="IDX1958"></A>
Non-rounded Akima spline with natural boundary conditions.  This method
uses the non-rounded corner algorithm of Wodicka.
</DL>

</P>
<P>
<DL>
<DT><U>Interpolation Type:</U> <B>gsl_interp_akima_periodic</B>
<DD><A NAME="IDX1959"></A>
Non-rounded Akima spline with periodic boundary conditions.  This method
uses the non-rounded corner algorithm of Wodicka.
</DL>

</P>
<P>
The following related functions are available:

</P>
<P>
<DL>
<DT><U>Function:</U> const char * <B>gsl_interp_name</B> <I>(const gsl_interp * <VAR>interp</VAR>)</I>
<DD><A NAME="IDX1960"></A>
This function returns the name of the interpolation type used by <VAR>interp</VAR>.
For example,

</P>

<PRE>
printf ("interp uses '%s' interpolation.\n", 
        gsl_interp_name (interp));
</PRE>

<P>
would print something like,

</P>

<PRE>
interp uses 'cspline' interpolation.
</PRE>

</DL>

<P>
<DL>
<DT><U>Function:</U> unsigned int <B>gsl_interp_min_size</B> <I>(const gsl_interp * <VAR>interp</VAR>)</I>
<DD><A NAME="IDX1961"></A>
This function returns the minimum number of points required by the
interpolation type of <VAR>interp</VAR>.  For example, Akima spline interpolation
requires a minimum of 5 points.
</DL>

</P>


<H2><A NAME="SEC402" HREF="gsl-ref_toc.html#TOC402">Index Look-up and Acceleration</A></H2>

<P>
The state of searches can be stored in a <CODE>gsl_interp_accel</CODE> 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.

</P>
<P>
<DL>
<DT><U>Function:</U> size_t <B>gsl_interp_bsearch</B> <I>(const double <VAR>x_array</VAR>[], double <VAR>x</VAR>, size_t <VAR>index_lo</VAR>, size_t <VAR>index_hi</VAR>)</I>
<DD><A NAME="IDX1962"></A>
This function returns the index i of the array <VAR>x_array</VAR> such
that <CODE>x_array[i] &#60;= x &#60; x_array[i+1]</CODE>.  The index is searched for
in the range [<VAR>index_lo</VAR>,<VAR>index_hi</VAR>].
</DL>

</P>
<P>
<DL>
<DT><U>Function:</U> gsl_interp_accel * <B>gsl_interp_accel_alloc</B> <I>(void)</I>
<DD><A NAME="IDX1963"></A>
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.
</DL>

</P>
<P>
<DL>
<DT><U>Function:</U> size_t <B>gsl_interp_accel_find</B> <I>(gsl_interp_accel * <VAR>a</VAR>, const double <VAR>x_array</VAR>[], size_t <VAR>size</VAR>, double <VAR>x</VAR>)</I>
<DD><A NAME="IDX1964"></A>
This function performs a lookup action on the data array <VAR>x_array</VAR>
of size <VAR>size</VAR>, using the given accelerator <VAR>a</VAR>.  This is how
lookups are performed during evaluation of an interpolation.  The
function returns an index i such that <CODE>x_array[i] &#60;= x &#60;
x_array[i+1]</CODE>.
</DL>

</P>
<P>
<DL>
<DT><U>Function:</U> void <B>gsl_interp_accel_free</B> <I>(gsl_interp_accel* <VAR>acc</VAR>)</I>
<DD><A NAME="IDX1965"></A>
This function frees the accelerator object <VAR>acc</VAR>.
</DL>

</P>


<H2><A NAME="SEC403" HREF="gsl-ref_toc.html#TOC403">Evaluation of Interpolating Functions</A></H2>

<P>
<DL>
<DT><U>Function:</U> double <B>gsl_interp_eval</B> <I>(const gsl_interp * <VAR>interp</VAR>, const double <VAR>xa</VAR>[], const double <VAR>ya</VAR>[], double <VAR>x</VAR>, gsl_interp_accel * <VAR>acc</VAR>)</I>
<DD><A NAME="IDX1966"></A>
<DT><U>Function:</U> int <B>gsl_interp_eval_e</B> <I>(const gsl_interp * <VAR>interp</VAR>, const double <VAR>xa</VAR>[], const double <VAR>ya</VAR>[], double <VAR>x</VAR>, gsl_interp_accel * <VAR>acc</VAR>, double * <VAR>y</VAR>)</I>
<DD><A NAME="IDX1967"></A>
These functions return the interpolated value of <VAR>y</VAR> for a given
point <VAR>x</VAR>, using the interpolation object <VAR>interp</VAR>, data arrays
<VAR>xa</VAR> and <VAR>ya</VAR> and the accelerator <VAR>acc</VAR>.
</DL>

</P>
<P>
<DL>
<DT><U>Function:</U> double <B>gsl_interp_eval_deriv</B> <I>(const gsl_interp * <VAR>interp</VAR>, const double <VAR>xa</VAR>[], const double <VAR>ya</VAR>[], double <VAR>x</VAR>, gsl_interp_accel * <VAR>acc</VAR>)</I>
<DD><A NAME="IDX1968"></A>
<DT><U>Function:</U> int <B>gsl_interp_eval_deriv_e</B> <I>(const gsl_interp * <VAR>interp</VAR>, const double <VAR>xa</VAR>[], const double <VAR>ya</VAR>[], double <VAR>x</VAR>, gsl_interp_accel * <VAR>acc</VAR>, double * <VAR>d</VAR>)</I>
<DD><A NAME="IDX1969"></A>
These functions return the derivative <VAR>d</VAR> of an interpolated
function for a given point <VAR>x</VAR>, using the interpolation object
<VAR>interp</VAR>, data arrays <VAR>xa</VAR> and <VAR>ya</VAR> and the accelerator
<VAR>acc</VAR>. 
</DL>

</P>
<P>
<DL>
<DT><U>Function:</U> double <B>gsl_interp_eval_deriv2</B> <I>(const gsl_interp * <VAR>interp</VAR>, const double <VAR>xa</VAR>[], const double <VAR>ya</VAR>[], double <VAR>x</VAR>, gsl_interp_accel * <VAR>acc</VAR>)</I>
<DD><A NAME="IDX1970"></A>
<DT><U>Function:</U> int <B>gsl_interp_eval_deriv2_e</B> <I>(const gsl_interp * <VAR>interp</VAR>, const double <VAR>xa</VAR>[], const double <VAR>ya</VAR>[], double <VAR>x</VAR>, gsl_interp_accel * <VAR>acc</VAR>, double * <VAR>d2</VAR>)</I>
<DD><A NAME="IDX1971"></A>
These functions return the second derivative <VAR>d2</VAR> of an interpolated
function for a given point <VAR>x</VAR>, using the interpolation object
<VAR>interp</VAR>, data arrays <VAR>xa</VAR> and <VAR>ya</VAR> and the accelerator
<VAR>acc</VAR>. 
</DL>

</P>
<P>
<DL>
<DT><U>Function:</U> double <B>gsl_interp_eval_integ</B> <I>(const gsl_interp * <VAR>interp</VAR>, const double <VAR>xa</VAR>[], const double <VAR>ya</VAR>[], double <VAR>a</VAR>, double <VAR>b</VAR>, gsl_interp_accel * <VAR>acc</VAR>)</I>
<DD><A NAME="IDX1972"></A>
<DT><U>Function:</U> int <B>gsl_interp_eval_integ_e</B> <I>(const gsl_interp * <VAR>interp</VAR>, const double <VAR>xa</VAR>[], const double <VAR>ya</VAR>[], double <VAR>a</VAR>, double <VAR>b</VAR>, gsl_interp_accel * <VAR>acc</VAR>, double * <VAR>result</VAR>)</I>
<DD><A NAME="IDX1973"></A>
These functions return the numerical integral <VAR>result</VAR> of an
interpolated function over the range [<VAR>a</VAR>, <VAR>b</VAR>], using the
interpolation object <VAR>interp</VAR>, data arrays <VAR>xa</VAR> and <VAR>ya</VAR> and
the accelerator <VAR>acc</VAR>.
</DL>

</P>


<H2><A NAME="SEC404" HREF="gsl-ref_toc.html#TOC404">Higher-level Interface</A></H2>

<P>
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
<CODE>gsl_interp</CODE> functions but maintain a copy of this data in the
<CODE>gsl_spline</CODE> object.  This removes the need to pass both <VAR>xa</VAR>
and <VAR>ya</VAR> as arguments on each evaluation. These functions are
defined in the header file <TT>`gsl_spline.h'</TT>.

</P>
<P>
<DL>
<DT><U>Function:</U> gsl_spline * <B>gsl_spline_alloc</B> <I>(const gsl_interp_type * <VAR>T</VAR>, size_t <VAR>size</VAR>)</I>
<DD><A NAME="IDX1974"></A>
</DL>

</P>
<P>
<DL>
<DT><U>Function:</U> int <B>gsl_spline_init</B> <I>(gsl_spline * <VAR>spline</VAR>, const double <VAR>xa</VAR>[], const double <VAR>ya</VAR>[], size_t <VAR>size</VAR>)</I>
<DD><A NAME="IDX1975"></A>
</DL>

</P>
<P>
<DL>
<DT><U>Function:</U> void <B>gsl_spline_free</B> <I>(gsl_spline * <VAR>spline</VAR>)</I>
<DD><A NAME="IDX1976"></A>
</DL>

</P>
<P>
<DL>
<DT><U>Function:</U> double <B>gsl_spline_eval</B> <I>(const gsl_spline * <VAR>spline</VAR>, double <VAR>x</VAR>, gsl_interp_accel * <VAR>acc</VAR>)</I>
<DD><A NAME="IDX1977"></A>
<DT><U>Function:</U> int <B>gsl_spline_eval_e</B> <I>(const gsl_spline * <VAR>spline</VAR>, double <VAR>x</VAR>, gsl_interp_accel * <VAR>acc</VAR>, double * <VAR>y</VAR>)</I>
<DD><A NAME="IDX1978"></A>
</DL>

</P>
<P>
<DL>
<DT><U>Function:</U> double <B>gsl_spline_eval_deriv</B> <I>(const gsl_spline * <VAR>spline</VAR>, double <VAR>x</VAR>, gsl_interp_accel * <VAR>acc</VAR>)</I>
<DD><A NAME="IDX1979"></A>
<DT><U>Function:</U> int <B>gsl_spline_eval_deriv_e</B> <I>(const gsl_spline * <VAR>spline</VAR>, double <VAR>x</VAR>, gsl_interp_accel * <VAR>acc</VAR>, double * <VAR>d</VAR>)</I>
<DD><A NAME="IDX1980"></A>
</DL>

</P>
<P>
<DL>
<DT><U>Function:</U> double <B>gsl_spline_eval_deriv2</B> <I>(const gsl_spline * <VAR>spline</VAR>, double <VAR>x</VAR>, gsl_interp_accel * <VAR>acc</VAR>)</I>
<DD><A NAME="IDX1981"></A>
<DT><U>Function:</U> int <B>gsl_spline_eval_deriv2_e</B> <I>(const gsl_spline * <VAR>spline</VAR>, double <VAR>x</VAR>, gsl_interp_accel * <VAR>acc</VAR>, double * <VAR>d2</VAR>)</I>
<DD><A NAME="IDX1982"></A>
</DL>

</P>
<P>
<DL>
<DT><U>Function:</U> double <B>gsl_spline_eval_integ</B> <I>(const gsl_spline * <VAR>spline</VAR>, double <VAR>a</VAR>, double <VAR>b</VAR>, gsl_interp_accel * <VAR>acc</VAR>)</I>
<DD><A NAME="IDX1983"></A>
<DT><U>Function:</U> int <B>gsl_spline_eval_integ_e</B> <I>(const gsl_spline * <VAR>spline</VAR>, double <VAR>a</VAR>, double <VAR>b</VAR>, gsl_interp_accel * <VAR>acc</VAR>, double * <VAR>result</VAR>)</I>
<DD><A NAME="IDX1984"></A>
</DL>

</P>


<H2><A NAME="SEC405" HREF="gsl-ref_toc.html#TOC405">Examples</A></H2>

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

</P>

<PRE>
#include &#60;stdlib.h&#62;
#include &#60;stdio.h&#62;
#include &#60;math.h&#62;
#include &#60;gsl/gsl_errno.h&#62;
#include &#60;gsl/gsl_spline.h&#62;

int
main (void)
{
  int i;
  double xi, yi, x[10], y[10];

  printf ("#m=0,S=2\n");

  for (i = 0; i &#60; 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 &#60; 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;
}
</PRE>

<P>
The output is designed to be used with the GNU plotutils
<CODE>graph</CODE> program,

</P>

<PRE>
$ ./a.out &#62; interp.dat
$ graph -T ps &#60; interp.dat &#62; interp.ps
</PRE>

<P>
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</CODE>.

</P>



<H2><A NAME="SEC406" HREF="gsl-ref_toc.html#TOC406">References and Further Reading</A></H2>

<P>
Descriptions of the interpolation algorithms and further references can
be found in the following book,

</P>

<UL>
<LI>C.W. Ueberhuber,

<CITE>Numerical Computation (Volume 1), Chapter 9 "Interpolation"</CITE>,
Springer (1997), ISBN 3-540-62058-3.
</UL>

<P>

</P>
<P><HR><P>
<p>Go to the <A HREF="gsl-ref_1.html">first</A>, <A HREF="gsl-ref_25.html">previous</A>, <A HREF="gsl-ref_27.html">next</A>, <A HREF="gsl-ref_50.html">last</A> section, <A HREF="gsl-ref_toc.html">table of contents</A>.
</BODY>
</HTML>