Алгоритм реализован в виде программы на языке С, которую вы можете скачать в разделе программы.

Опишем сначала, что такое интерполяция и как эта задача решена в данном случае.

Пусть на некотором отрезке в точках x₀, x₁, x₂, ... x_N нам известны значения некоторой функции f(x), а именно y₀, y₁, y₂, ... y_N.

Требуется построить интерполирующую функцию F(x), такую, что она принимает в указаных точках те же значения, т.е. F(x₀) = y₀, F(x₁) = y₁, ... F(x_N) = y_N.

Геометрически это значит, что нужно найти кривую y = F(x) определенного типа, проходящую через систему заданных точек. В такой общей постановке задача может иметь бесчисленное множество решений или совсем не иметь решений. Однако эта задача становится однозначной, если вместо произвольной функции F(x) искать полином P_N(x), степени N, удовлетворяющий описаным выше условиям, т.е. который проходит через все заданные точки.

Если искать функцию в виде F(x) = P_N(x) = a₀+ a₁x¹+ a₂x²+ ... a_Nx^N то эта интерполяция называется интерполяцией каноническим полиномом.

Выбор многочлена степени N основан на том факте, что через N+1 точку проходит единственная кривая степени N. Подставив значения функции y_i в точках x_i в выражение полинома получим систему N+1 линейных уравнений.

a0 +	a1x0 +	a2x02 +	...	aNx0N = y0
a0 +	a1x1 +	a2x12 +	...	aNx1N = y1
a0 +	a1x2 +	a2x22 +	...	aNx2N = y2
...
a0 +	a1xN +	a2xN2 +	...	aNxNN = yN

В данной системе известны все y_i, x_i^k, неизвестны все a_i. Число неизвестных равно числу уравнений, а именно N+1. Уравнения являются линейными относительно a_i, следовательно их можно решить, а решение существует и единственно (при некоторых условиях на детерминант системы уравнений).

Для решения системы воспользуемся методом Гаусса в котором система уравнений приводится к треугольному виду. Я уже описывал этот метод. Описываемая в данном разделе программа основана на уже упомянутой, поэтому я не буду подробно останавливаться на сущности метода Гаусса, лишь приведу текст программы и кратко опишу методику работы с ней, а также ее структуру, а именно, для чего служит такое большое количество подпрограмм в этой программе.

//Classic polinomial interpolation
//when we have two columns of data x and y in input file:
//
//x0 y0
//x1 y1
//...
//xn yn
//
//and we want to find such function f(x) = a0 + a1*x + a2*x^2 + ... + an*x^n
//where f(xi) = yi
//
//the result is ai column
#include <stdio.h>
#include <process.h>
#include <math.h>
float *a, **x, *y;
int N;
char filename[256];
FILE* InFile=NULL;
void count_num_lines(){
   //count number of lines in input file - number of equations
   int nelf=0;       //non empty line flag
   do{
       nelf = 0;
       while(fgetc(InFile)!='\n' && !feof(InFile)) nelf=1;
       if(nelf) N++;
   }while(!feof(InFile));
}

void freematrix(){
   //free memory for matrixes
   int i;
   for(i=0; i<N; i++){
       delete [] x[i];
   }
   delete [] a;
   delete [] x;
   delete [] y;
}

void allocmatrix(){
   //allocate memory for matrixes
   int i,j;
   a = new float[N];
   y = new float[N];
   x = new float*[N];
   if(x==NULL || y==NULL || a==NULL){
       printf("\nNot enough memory to allocate for %d equations.\n", N);
       exit(-1);
   }
   for(i=0; i<N; i++){
       x[i] = new float[N];
       if(x[i]==NULL){
	   printf("\nNot enough memory to allocate for %d equations.\n", N);
       }
   }
   for(i=0; i<N; i++){
       for(j=0; j<N; j++){
	   x[i][j]=0;
       }
       y[i]=0;
       a[i]=0;
   }
}

void readmatrix(){
   int i=0,j=0;
   //read matrixes from input file
   for(i=0; i<N; i++){
       fscanf(InFile, "%f", &x[i][1]);
       fscanf(InFile, "%f", &y[i]);
   }
   //init square x matrix
   for(i=0; i<N; i++){
       for(j=0; j<N; j++){
           x[i][j] = pow(x[i][1], j);
       }
   }
}

void printmatrix(){
   //print system of equations
   int i=0,j=0;
   printf("\n");
   for(i=0; i<N; i++){
       for(j=0; j<N; j++){
	   printf(" %+f*A%d", x[i][j], j);
       }
       printf(" =%f\n", y[i]);
   }
}

void testsolve(){
   //test that ax=y
   int i=0, j=0;
   printf("\n");
   for(i=0; i<N; i++){
       float s = 0;
       for(j=0; j<N; j++){
	   s += x[i][j]*a[j];
       }
       printf("%f\t%f\n", s, y[i]);
   }
}
void printresult(){
   int i=0;
   printf("\n");
   printf("Result\n");
   for(i=0; i<N; i++){
       printf("A%d = %f\n", i, a[i]);
   }
}
void diagonal(){
   int i, j, k;
   float temp=0;
   for(i=0; i<N; i++){
       if(x[i][i]==0){
	   for(j=0; j<N; j++){
	       if(j==i) continue;
	       if(x[j][i] !=0 && x[i][j]!=0){
		   for(k=0; k<N; k++){
		       temp = x[j][k];
		       x[j][k] = x[i][k];
		       x[i][k] = temp;
		   }
		   temp = y[j];
		   y[j] = y[i];
		   y[i] = temp;
		   break;
	       }
	   }
       }
   }
}
void cls(){
   for(int i=0; i<25; i++) printf("\n");
}
void main(){
   int i=0,j=0, k=0;
   cls();
   do{
       printf("\nInput filename: ");
       scanf("%s", filename);
       InFile = fopen(filename, "rt");
   }while(InFile==NULL);
   count_num_lines();
   allocmatrix();
   rewind(InFile);
   //read data from file
   readmatrix();
   //check if there are 0 on main diagonal and exchange rows in that case
   diagonal();
   fclose(InFile);
   printmatrix();
   //process rows
   for(k=0; k<N; k++){
       for(i=k+1; i<N; i++){
	   if(x[k][k]==0){
	       printf("\nSolution is not exist.\n");
	       return;
	   }
	   float M = x[i][k] / x[k][k];
	   for(j=k; j<N; j++){
	       x[i][j] -= M * x[k][j];
	   }
	   y[i] -= M*y[k];
       }
   }
   printmatrix();
   for(i=N-1; i>=0; i--){
       float s = 0;
       for(j = i; j<N; j++){
	   s = s+x[i][j]*a[j];
       }
       a[i] = (y[i] - s) / x[i][i];
   }

   InFile = fopen(filename, "rt");
   readmatrix();
   fclose(InFile);
   printmatrix();
   testsolve();
   printresult();
   freematrix();
}

Чтобы воспользоваться этой программой, нужно запустить скомпилированный исполняемый файл. В первую очередь программа спросит, откуда ей брать данные для интерполяции. Создайте в любом текстовом редакторе (но только не в Word-е а, например в notepad-е) файл, где напишите значения x_i, y_i, построчно через пробел, приблизительно так:

x0	y0
x1	y1
x2	y2
...
xN	yN

Например:

1	5
2	3
3	2.5
4	2
5	0

Этот файл необходимо создать в той директории, где лежит программа, иначе она его не найдет. В результате работы программы, она выдаст нечто вроде:

Result
A0 = 10.000000
A1 = -7.000000
A2 = 2.250000
A3 = -0.250000
A4 = 0.000000

Это и есть решение системы уравнений, т.е. набор неизвестных коэффициентов a_i. Т.е. по нашим данным был построен классический интерполяционный многочлен

F(x) = 10 - 7x + 2.25x² - 0.25x³.

Так случайно получилось, что он третьей степени, хотя для пяти точек в общем виде должен получаться многочлен четвертой степени.

Теперь для того, чтобы узнать значение функции в некоторой промежуточной точке, например x=2.5, которой не было в таблице, достаточно подставить значение x=2.5 в формулу:

F(2.5) = 10 - 7*2.5 + 2.25*2.5² - 0.25*2.5³ = 2.65525

что является правдоподобной оценкой значения функции.

Коротко опишем, для чего служит такое большое количество подпрограмм в данной программе:

• void count_num_lines() - подсчитывает количество точек, где задана функция
• void allocmatrix() - выделяет память для массивов a, x, y
• void readmatrix() - прочитывает из файла координаты точек и вычисляет x_i^k
• void printmatrix() - распечатывает систему уравнений
• void diagonal() - делает так, чтобы на главной диагонали не было нулей, чтобы не пришлось однажды делить на ноль в процессе приведения системы к треугольному виду
• void testsolve() - подставляет решение в систему и распечатывает рядом то, что получается в левой части уравнения и для сравнения печатает рядом свободные члены, те, что в правой части уравнения; получившиеся два столбца должны совпадать, если уравнения решены правильно
• void printresult() - распечатывает получившийся столбец решений
• void freematrix() - освобождает память, которая была выделена ранее
• void cls() - стирает экран в начале работы программы
• void main() - основная функция из которой последовательно вызываются все вышеперечисленные функции, и проходит процесс приведения системы уравнений к треугольному виду и обратная прогонка.

Надеюсь, что вам стало все понятно.