Bessel’s Interpolation in C++

Interpolation is a type of estimation technique of unknown value which lies between know values. Interpolation is the process of constructing new data points between the range of a discrete set of know data points.

An application or reason to use interpolation is that it might reduce computation costs. When the formula (function) to calculate certain values is too complicated or costly to compute, we prefer using interpolation. A few data points are calculated using the original function, the rest of them can be estimated using interpolation. These may not be completely accurate but fairly close!

So basically here the reduced computation cost and simplicity is outweighing the loss from interpolation error.

Bessel’s Interpolation Formula

f(u) = {(f(0)+f(1))/2} + {u - ½}𝛥f(0) + {u(u-1)/2!}{(𝛥2 f(-1) + 𝛥2 f(0))/2} + {u(u-1)(u - ½)/3!}𝛥3f(-1) + {u(u+1)(u-1)(u-2)/4!}{(𝛥4f(-2) + 𝛥4f(-1))/2}+..

here,

f(0) is origin point which is usually the mid point.

u = x - f(0) / h, gh is interval of difference

Example

Program to illustrate Bassel’s interpolation −

 Live Demo

#include <iostream>
using namespace std;
float calU(float u, int n){
   if (n == 0)
      return 1;
   float result = u;
   for (int i = 1; i <= n / 2; i++)
      result = result*(u - i);
   for (int i = 1; i < n / 2; i++)
      result = result*(u + i);
   return result;
}
int factorial(int n){
   if(n == 1)
      return 1;
   return n * factorial(n-1);
}
int main(){
   int n = 6;
   float x[] = { 50, 51, 52, 53, 54, 55 };
   float y[n][n];
   y[0][0] = 8.000;
   y[1][0] = 7.746;
   y[2][0] = 7.674;
   y[3][0] = 7.571;
   y[4][0] = 7.469;
   y[5][0] = 7.231;
   for (int i = 1; i < n; i++)
      for (int j = 0; j < n - i; j++)
         y[j][i] = y[j + 1][i - 1] - y[j][i - 1];
   float value = 53.2;
   float sum = (y[2][0] + y[3][0]) / 2;
   int index;
   if (n % 2)
      index = n/2;
   else
      index = n/2 - 1;
   float u = (value - x[index]) / (x[1] - x[0]);
   for (int i = 1; i < n; i++) {
      if (i % 2)
         sum+= (((u-(0.5))*calU(u, i - 1)*y[index][i])/factorial(i));
      else
         sum+= ((calU(u, i)*(y[index][i]+y[-- index][i])/(factorial(i)*2)));
   }
   cout<<"Value at "<<value<<" found using Bessels's interpolation is "<<sum;
   return 0;
}

Output

Value at 53.2 found using Bessels's interpolation is 7.54985