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: 
    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. No range checking is done either.

    You need to include <stdio.h> and <stdlib.h> for printf, malloc, and rand

  2. gsl_vector: 
    #define RND (double)rand()/RAND_MAX
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 know their sizes and by default do range-checking. Rrange-checking can be disabled by recompiling your program with the preprocessor definition GSL_RANGE_CHECK_OFF.

    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`

  3. Your own vector library

    Something like here.

C++
  • std::vector<double>

    This is a useful container from the Standard Template Library (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;
for(auto &e : v) cout << e <<endl; /* C++11 */
delete w;

    (you might need -std=c++11 to compile the range-based for-loop).
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]);
for(double e : v) System.out.println(e); /* 1.5.0 */

However, your own vector/matrix library (like this) may still come usefull as you can define some useful methods like add, subtract, scale, norm etc.

Python

Python has lists as containers although for numerical computaions one often uses arrays from NumPy. However, since an interpreted language is going to be slow anyway you can just use the built-in lists. 

n=5
v=[ i for i in range(n-1) ]
for e in v: print(e)
for i, e in enumerate(v): print(i,e)

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. For example, for C-arrays: 
      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 function, typically called qspline_alloc which allocates the necessary memory, calculates the spline coefficients b and c and returns the (pointer to the) properly filled qspline structure: 

      struct qspline* qspline_alloc(int n, double* x,double* y);
      Tthe name of this function usually ends with _alloc to underline that it allocates memory in the heap which then needs to be freed.

    3. Implement a function 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. Install it with sudo apt-get install plotutils. An example of graph usage can be found here and here. You can also use 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 (integral) of a tabulated function by building a quadratic spline and then analytically differentiating and integrating the spline.

    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

    ∫x0z S(x)dx = ∑k<i ∫xkxk+1 Sk(x)dx + ∫xiz Si(x)dx , where x[i]≤z≤x[i+1]

  3. (1 points) Cubic spline

    • Implement cubic spline with derivatives and integrals.
    • Make the first/second derivatives at the end-points input parameters to cubic spline.

  4. (0 points) Comparison with library routines

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