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

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

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

Требуется определить параметры a_i многочлена вида
F(x) = a₀ + ax + a₂x² + ... + a_kx^k, где k<N
такого, что сумма квадратов отклонений значений y от значений функции F(y) в заданных точках x была минимальной, т.е.
S=Σ [y_i - F(x_i, a₀, a₁...a_k)]² -> min

Геометрически это значит, что нужно найти кривую y = F(x), полином, который проходит как можно ближе к каждой из заданной точек.

Такая задача может быть решена, если решить систему уравнений вида:

a0n +	a1Σ xi +	a2Σ xi2 +	... +	akΣ xik =	Σ yi
a0Σ xi +	a1Σ xi2 +	a2Σ xi3 +	... +	akΣ xik+1 =	Σ xiyi
...
a0Σ xik +	a1Σ xik+1 +	a2Σ xik+2 +	... +	akΣ xi2k =	Σ xikyi

где везде под символом Σ подразумевается суммирование по i от 0 до N.

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

#include <stdio.h>
#include <process.h>
#include <math.h>
float *a, *b, *x, *y, **sums;
int N, K;
//N - number of data points
//K - polinom power
//K<=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<K+1; i++){
       delete [] sums[i];
   }
   delete [] a;
   delete [] b;
   delete [] x;
   delete [] y;
   delete [] sums;
}

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

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

void printresult(){
   //print polynom parameters
   int i=0;
   printf("\n");
   for(i=0; i<K+1; i++){
       printf("a[%d] = %f\n", i, a[i]);
   }
}
void diagonal(){
   int i, j, k;
   float temp=0;
   for(i=0; i<K+1; i++){
       if(sums[i][i]==0){
	   for(j=0; j<K+1; j++){
	       if(j==i) continue;
	       if(sums[j][i] !=0 && sums[i][j]!=0){
		   for(k=0; k<K+1; k++){
		       temp = sums[j][k];
		       sums[j][k] = sums[i][k];
		       sums[i][k] = temp;
		   }
		   temp = b[j];
		   b[j] = b[i];
		   b[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();
   printf("\nNumber of points: N=%d", N);
   do{
       printf("\nInput power of approximation polinom K<N: ");
       scanf("%d", &K);
   }while(K>=N);
   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<K+1; k++){
       for(i=k+1; i<K+1; i++){
	   if(sums[k][k]==0){
	       printf("\nSolution is not exist.\n");
	       return;
	   }
	   float M = sums[i][k] / sums[k][k];
	   for(j=k; j<K+1; j++){
	       sums[i][j] -= M * sums[k][j];
	   }
	   b[i] -= M*b[k];
       }
   }
   //printmatrix();
   for(i=(K+1)-1; i>=0; i--){
       float s = 0;
       for(j = i; j<K+1; j++){
	   s = s + sums[i][j]*a[j];
       }
       a[i] = (b[i] - s) / sums[i][i];
   }
   printresult();
   freematrix();
}

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

x0	y0
x1	y1
x2	y2
...
xN	yN

Например:

1	5.95
2	20.95
3	51.9
4	105
5	186
6	301
7	456.1
8	657.1

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

Программа спросит степень полинома, каким вы хотите аппроксимировать данные. При этом степень полинома должна быть менее числа заданных точек (в данном случае восьми). Введите, например, 3. В результате работы программы, она выдаст нечто вроде:

a[0] = 0.996902
a[1] = 1.938750
a[2] = 2.016942
a[3] = 0.999054

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

F(x) = 1 + 2x + 2x² + x³. (я округлил решение до целых чисел)

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

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

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