Interpolation

Here we shall go through a simple example on interpolation. We wish to interpolate to find the value of the function $y(x)$ (which has no formula) at an arbitrary position $a$ given the following data:

i	$x_i$	$y(x_i)$
0	0	1.90246
1	0.4	1.50246
2	0.8	1.10507
3	1.2	0.739343
4	1.6	0.454018
5	2	0.262869
6	2.4	0.146872
7	2.8	0.0799558
8	3.2	0.0414204
9	3.6	0.0175786
10	4	0

• The first task is to write a program and input the data here into two vectors $x$ and $y$.
• Then we wish to write a function interpolate that will take the two vectors, along with the value $a$ and an integer to specify precision as input, then return the interpolated value $y(x=a)$.
• The method we shall use to interpolate is a simple method using a Lagrange polynomial.

Prototype for the function

The function prototype should look like:

double interpolate(vector< double>& x,vector< double>& y, double a, int degree)
{
   // interpolate in here...
}

where $x$ and $y$ are vectors containing the data, $a$ is the position at we wish to evaluate the function, and degree is the number of points that we wish to use.

Generating a Lagrange polynomial

The Lagrange polynomial is defined as the linear combination:

$$L(a) = \sum_{j=0}^{n} y_jl_j(a)$$

where the polynomial $l$ is defined by:

$$l_j(a) = \Pi_{i=0,i\neq j}^{n} \frac{a - x_i}{x_{j} - x_i}$$

Finding where to interpolate

The function is to be split into two parts:

• finding the value of $x_i$ closest to $a$,
• interpolating using the closest points to $a$

The idea here is that we can let degree be smaller than the size of the vectors $x$ and $y$.

In order to find which points to use in the interpolation, let us assume that (as is the case above, and when interpolating finite difference grids)

$$x_i = i\Delta x.$$

The value of $\Delta x$ here is just the difference between any two nodes in $x$. Then given the value $a$, we know that there exists some $i^*$ such that

$$a \geq i^*\Delta x,\quad \mathrm{and} \quad a < (i^*+1)\Delta x$$

To find $i^*$ with using C++ we can simply say

int istar;
istar = a/dx;


However, if we wish to find the closest node to $a$ then write

int istar;
istar = int(a/dx+0.5);


Try this out to see what the subtle difference between the two is.

Example - Linear Interpolation

Given the protoype for the function as given above, we shall choose degree equal to 2 as an example. There are two stages:

• (1) Given the value $a$, find $i^*$ (use the method to find the node below)
  
  int istar;
  istar = a/dx;
  
  
  
• (2) Evaluate the Lagrange polynomial at $x=a$
  
  double sum=0.;
  sum = sum + (a - x[istar+1])/(x[istar]-x[istar+1])*y[istar];
  sum = sum + (a - x[istar])/(x[istar+1]-x[istar])*y[istar+1];
  return sum;
  
  
  
• (3) Test and debug the code!!! You can check that the polynomial above is the equation of a line passing throught the points $(x_{i^*},y_{i^*})$ and $(x_{i^*+1},y_{i^*+1})$.

Higher degree approximations

• After testing the code, see if you can write the code into loops (see wiki).
• Depending on whether degree is odd or even, it will be more appropriate to choose the nearest node or the node below $a$. Think about which one is needed in each case and why.
• Write the function for the general case and test it with degree equal to 4. The results above are from the European put example from the previous tutorial on Crank-Nicolson. Set $a=1.8$ and see how accurately the results agree with the Black-Scholes formula.
• Be careful! Don't choose degree too high. It should never need to be greater than 5.