Homework "Splines"

Objective

Hints

You can use either our "vector" class or "double[]" to keep the x and y vectors.

The exercise (parts B and C) can be solved using different programming styles:

1. Procedural programming: One makes two functions (procedures). The first one, say

   static (vector,vector) build_qspline(vector x,vector y){
   	/* build and return (b,c); */
   	}
   

   builds the (quadratic) spline-coefficients b and c. The second one, say

   static double evaluate_qspline(vector x,vector y,vector b,vector c,double z){
   	/* evaluate si(z) using x,y,b,c and return */
   	}
   

   evaluates the spline at point "z" using the vectors x, y, b, c. The disadvantage is that one has to explicitly allocate and pass around the spline-coefficients b and c:

   (vector b,vector c) = build_qspline(x,y);
   double fz = evaluate_qspline(x,y,b,c,z);
   

2. Object oriented programming: One makes an object that holds both the data and the appropriate functions (now called member functions or methods). The methods have direct access to the data in the structure, for example,

   public class qspline{
   	vector x,y,b,c;
   	public qspline(vector xs, vector ys){
   		/* copy xs,ys into x,y and build b,c */
   		}
   	public double evaluate(double z){
   		/* evaluate the spline using x,y,b,c */
   		}
   }
   

   Now the procedural style is replaced with the object-oriented style,

   var myspline = new qspline(x,y);
   double fz = myspline.evaluate(z);
   

   No need to pass the vectors x,y,b,c to "myspline.evaluate" – they are kept in the "myspline" instance of the class.

3. Functional programming: The function for spline-evaluation is created and returned. This returned function captures all the needed information in itself (in its environment, to be precise). For example,

   static Func<double,double> qspline(vector xs,vector ys){
   	vector x, y, b, c; /* will be captured when qspline returns */
   	/* copy xs,ys into x,y and calculate b and c */
   	return delegate(double z){/* evaluate and return spline(z) */};
   }
   

   Now the procedural style is replaced by the functional style,

   Func<double,double> myspline = qspline(x,y);
   double fz=myspline(z);
   

   No need to pass the vectors x,y,b,c to "myspline" – they are captured in the (environment of the) function "myspline".

Tasks

1. (6 points) Linear spline (linear interpolation)
   1. Implement a function that makes linear spline interpolation from a table {x[i],y[i]} at a given point z. Here you don't need to either pre-calculate anything or use any structs. The signature of your function should be something like

      public static double linterp(double[] x, double[] y, double z){
              int i=binsearch(x,z);
              double dx=x[i+1]-x[i]; if(!(dx>0)) throw new Exception("uups...");
              double dy=y[i+1]-y[i];
              return y[i]+dy/dx*(z-x[i]);
              }
      

      Location of the index i of the interval containing z (such that x[i]<z<x[i+1]), must be done using binary search, like this:

      public static int binsearch(double[] x, double z)
      	{/* locates the interval for z by bisection */ 
      	if( z<x[0] || z>x[x.Length-1] ) throw new Exception("binsearch: bad z");
      	int i=0, j=x.Length-1;
      	while(j-i>1){
      		int mid=(i+j)/2;
      		if(z>x[mid]) i=mid; else j=mid;
      		}
      	return i;
      	}
      

   2. Implement a function that calculates the integral of the linear spline from the point x[0] to the given point z. The integral must be calculated analytically as it is the integral of a linear function. The signature could be something like

      double linterpInteg(double[] x, double[] y, double z)
      

   3. Make some indicative plots to prove that your linear spline and your integrator work as intended. For example, take the table

      {xi=0,1,…,9; yi=cos(xi)},
      

      and plot its linear interpolant together with interpolant's anti-derivative.

2. (3 points) Quadratic spline

   Implement quadratic spline with derivative and definite integral (anti-derivative). Note that quadratic spline is only intended for learning, for practical applications always use the cubic spline.

   In C-sharp a quadratic spline can be kept in a class like the following,

   public class qspline {
   	vector x,y,b,c;
   	public qspline(vector xs,vector ys){
   		/* x=xs.copy(); y=ys.copy(); calculate b and c */
   		}
   	public double evaluate(double z){/* evaluate the spline */}
   	public double derivative(double z){/* evaluate the derivative */}
   	public double integral(double z){/* evaluate the integral */}
   	}
   

   Alternatively, you can use functional programming style:

   static Func<double, int, double> qspline(vector xs,vector ys){
   	/* x=xs.copy(); y=ys.copy(); calculate b and c */
   	return delegate(double z,int deriv=0){
   		if(deriv==1) /* return derivative */
   		if(deriv==-1) /* return integral */
   		else /* return spline */
   		}; /* x,y,b,c are captured here */
   }
   

   Make some indicative plots to prove that your linear spline and your integrator work as intended. For example, {x_i=0,1,…,9; y_i=sin(x_i)}.

   Hint: you can debug your quadratic-spline by considering the following {x_i,y_i} tables,
   {x_i=i, y_i=1} , i=1…5
   {x_i=i, y_i=x_i} , i=1…5
   {x_i=i, y_i=x_i²} , i=1…5
   Calculate manually the parameters {b_i, c_i} of the corresponding quadratic-splines, and compare the results with your quadratic-spline program.

3. (1 points) Cubic spline

   Implement the cubic spline with derivative and definite integral (anti-derivative).

   Check that the built-in cubic splines in pyxplot/gnuplot (or the spline utility from plotutils) produces a similar cubic spline to your implementation.