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)

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

Programming styles

The exercise can be solved using different styles of programming:

  1. Procedural programming (typically C, Pascal, Fortran):

    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. Implement 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. Implement 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);     /* C or procedural C++ */
function<double(double)> lspline(vector&x,vector&y);       /* functional C++ */
def linterp(x:list, y:list, z:float):                      # procedural python
    2. Implement the following functions/objects/methods for quadratic spline interpolation,
      • C : 
        struct qspline * qspline_alloc(int n, double *x, double *y); /* builds the quadratic spline */
double qspline_evaluate(struct qspline *s, double z);        /* evaluates the spline at point z */
double qspline_derivative(struct qspline *s, double z); /* evaluates the derivative of the spline at point z */
double qspline_integral(struct qspline *s, double z);  /* evaluates the integral of the spline from x[0] to z */
void qspline_free(struct qspline *s); /* free memory allocated in qspline_alloc */
      • C++ : 
        function<double(double,int)> qspline(vector<double>& x,vector<double>& y); /* functional style */
        or 
        struct qspline { /* object-oriented style */
   vector<double> x,y,b,c;
   qspline(vector<double> xtable, vector<double> ytable);
   double eval(double z);
   double deriv(double z);
   double integ(double z);
};
      • Python : 
        def qspline(x:list,y:list) :
	# calculate b[i],c[i]
	def eval(z:float, deriv:int=0)
		if deriv==1 : # calculate and return derivative
		elif deriv==-1 : # calculate and return integral
		else : # calculate and return spline
	return eval # functional style

    Hints:

    1. To illustrate that your routines work as intended, tabulate some known function and interpolate the table using your routines, including derivative and integral. Compare with the exact results.
    2. Location of the index i: x[i]<z<x[i+1] must(!) be done using the binary search.
    3. 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 if you like.

  2. (3 points) Cubic spline

    • Implement cubic spline with derivatives and integrals.

  3. (1 points) Vector structure/class/object

    • Build you own structure/class with the corresponding functions/methods to work with vectors. In the minimum you have to implement the allocator/constructor; deallocator/destructor (unless you use a garbage collector); and getter/setter.
      • C : 
        typedef struct {int size; double* data;} vector;
vector* vector_alloc(int n); 
void    vector_free(vector* v);
double  vector_get(vector* v, int i);           /* v_i */
void    vector_set(vector* v, int i, double x); /* v_i = x */
      • C++ : 
        namespace num{
struct vec{
	int size;
	std::vector<double> data;
	vec(int n) : size(n), data(n) {};
	double get(int n)             {return data[i];};
	void set(int i, double value) {data[i]=value;};
	double& operator[] (int i)    {return data[i];};
	double& operator() (int i)    {return data[i];};
}
}
      • Python : 
        import math, array
class vector(object):
	def __init__(self, n:int, t:type='d') :
		self.size=n
		self.data=array.array(t,[0.0]*n)

	def get(self,i:int) : return self.data[i]
	def set(self,i:int,value): self.data[i]=value

	def __getitem__ (self,i:int) : return self.data[i]
	def __setitem__ (self,i:int,value) : self.data[i]=value
    • An example of the vector class in Python can be found here.
    • Add some useful functions/method to your structure/class and/or overload some useful operators.

  4. (0 points) Extra

    • Check that the spline utility from plotutils produces a similar cubic spline to your implementation.
    • For C,C++: check that the GSL cubic spline functions 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.
    • Make the first/second derivatives at the end-points input parameters to cubic spline.