Exercise "Interpolation"

Objective:

Implement functions for interpolation of tabulated data.

Introduction:

Arrays

A mathematical set {xi | i=1,...,n} can be represented in programming languages in many ways, usually (but not necessarily) as arrays with syntax x[i] for the element xi.

For example:

C, C++
  1. C-array: 
    #define RND (double)rand()/RAND_MAX
int n=5;
double v[n]; /* in the stack */
double * w = malloc(n*sizeof(double)); /* in the heap */
for(int i=0;i<n;i++) { v[i]=i, w[i]=v[i]; }
for(int i=0;i<n;i++) printf("%g %g\n",v[i],w[i]);
free(w); /* there should be a free for every alloc */

    Note that these arrays are passed as pointers (of the type double*) and do not remember their sizes. The size of the array must be passed together with the array.

    (You need to include "stdio.h" and "stdlib.h" for printf, malloc, and rand)

  2. gsl_vector: 
    int n=5;
gsl_vector * x = gsl_vector_alloc(n);
for(int i=0; i<x->size; i++) gsl_vector_set(x,i,RND);
for(int i=0; i<x->size; i++) printf("%g\n",gsl_vector_get(x,i));
gsl_vector_free(x);

    GSL vectors remember their sizes. You only need to pass the pointer to the vector.

    You need to i) install GSL (sudo apt-get install libgsl0-dev); ii) #include<gsl/gsl_vector.h>, and iii) put something like this into your makefile 

    CFLAGS = -Wall -std=gnu99 `gsl-config --cflags`
CXXFLAGS = -Wall -std=c++0x `gsl-config --cflags`
LDLIBS = -lstdc++ `gsl-config --libs`

C++
  • std::vector<double>

    This is a useful container from STL: 

    int n=5;
vector<double> v(n); /* in the stack */
auto w = new vector<double>(n); /* in the heap */
for(uint i=0; i<v.size(); i++) v[i]=RND;
for(uint i=0; i<v.size(); i++) (*w)[i]=v[i];
for(uint i=0; i<(*w).size(); i++) cout << (*w)[i] <<endl;
delete w;

Java

In Java arrays remember their sizes: 

int n=5;
double[] v= new double[n];
for(int i=0; i<v.length; i++) v[i]=Math.random();
for(int i=0; i<v.length; i++) System.out.println(v[i]);

Programming styles

The exercise can be solved using different styles of programming:

  1. Procedural programming (C,Pascal,Fortran-77,...):

    1. Define a structure which holds all the needed information about the quadratic spline: 
      struct qspline {int n; double *x,*y,*b,*c};
      • Hint: alternatively, you can define a new type qspline 
        typedef struct {int n; double *x, *y, *b, *c;} qspline;
        and then in the following instead of struct qspline simply use qspline.

    2. Implement a routine which allocates the necessary memory, calculates the spline coefficients c and returns the (pointer to the) properly filled qspline structure (the name of this function usually ends with _alloc to underline that it allocates memory in the heap which then needs to be freed): 

      struct qspline* qspline_alloc(int n, double* x,double* y);

    3. Implement a routine which—given the precalculated qspline structure—computes the spline interpolation function at a given point z: 

      double qspline_evaluate(struct qslpine* s, double z);
    4. Define a function, qspline_free, which frees the memory occupied by structure qspline. 
      void qspline_free(struct qslpine* s);

  2. Object oriented programming (C++, Java, Python, C#, ...): the structure (or class) holds not only data but also all the appropriate functions (called member functions) which have direct access to the data in the structure:

    1. Define a structure (class) qspline which holds all the data for spline; the construction function becomes the constructor method; the evaluation function becomes the member function; the qspline_free function becomes the destructor method which is defined automatically (I think). In C++: 
      struct qspline {
vector<double> x,y,b,c; /* the data for the spline */
qspline(vector<double> x, vector<double> y); /* constructor */
double evaluate(double z); /* evaluates the interpolation function */
};
    The std::vector<double> is a useful container (which remembers its length and can do bounds checking) in C++.

  3. Functional programming (C++11, Lisp, Haskell, Python, JavaScript, ...): only the function for spline-evaluation is created and returned. This function captures all the needed information in itself (in its environment, to be precise):

    1. Define a function which calculates the spline coefficients and returns a function which evaluates the spline interpolation function using the captured coefficients (and other data): 
      function<double(double)> construct_qspline_evaluation_function(vector<double> x, vector<double> y);
Tasks:
  1. (6 points) Linear and quadratic splines

    1. Implement a function which makes linear spline interpolation from a table {x[i], y[i]} at a given point z. Here you don't need to precalculate anything or use any structs. The signature of your function should be something like

      double linterp(int n, double *x, double *y, double z);

    2. Implement the following (self-evident) routines for quadratic spline interpolation, 
      struct qspline * qspline_alloc(int n, double *x, double *y);
double qspline_evaluate(struct qspline *s, double z);
void qspline_free(struct qspline *s);
    3. Perform linear spline and quadratic spline of the follolwing data (or, even better, your own illustrative data set),

      for(int i = 0; i < 10; i++) { x[i] = i + 0.5 * sin (i); y[i] = i + cos (i * i); }

      and plot the two interpolating functions (together with the data points).

    Hints:

    1. Location of the index i: x[i]<z<x[i+1] must(!) be done using the binary search.
    2. The plot can be done using the graph utility from plotutils (an example of its usage can be found here), or using pyxplot or gnuplot.

  2. (3 points) Derivative and anti-derivative of tabulated function using quadratic spline

    Implement a function which estimates the derivative and the anti-derivative (indefinite integral) of a tabulated function by building a quadratic spline and then differentiating and integrating the spline (analytically!).

    The signatures of the functions should be like 

    double qspline_derivative(struct qspline *s, double z) /* returns s'(z) */
double qspline_integral(struct qspline *s, double z) /* returns \int_{x[0]}^z s(x) dx */

    Hint: the integration function should return
    z
    ∫ S(x)dx
    x0

  3. (1 points) Piecewise polynomial interpolation

    Implement a function which makes a polynomial intepolation using not all the points in the table (that would be susceptible to the Runge's phenomenon) but rather only few neighboring points to the right and to the left of the point of interpolation.

    The signature of the function should be something like 

    double piecewise_polynomial(int n, double x[], double y[], double z, int number_of_neighbors_to_use)

    Check your implementation on some indicative data of your choice.

    Implement also the derivatives and anti-derivatives.

    Investigate (analytically) whether this piecewise polynomial has continuous derivatives.

  4. (0 points) Cubic spline and comparison with library routines
    • Implement the cubic spline with derivatives and anti-derivatives.
    • For C,C++: check that the GSL cubic spline functions give the same result as your cubic spline.
    • For C++: make a C++ wrapper for the GSL cubic spline functions and check that they give the same result as your cubic spline.
    • For languages with available libraries: check that library implementation of cubic spline gives the same result as your spline.
    • For languages without good libraries: check that Octave cubic spline routine gives the same result as your implementation. On molveno Octave is simple octave, on lifa Octave is /usr/users/fedorov/bin/octave.