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

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

Пусть на некотором отрезке в точках 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) состоит из кусочков, а именно, на каждом из отрезков [x_k-1; x_k] функция F(x) является кубическим полиномом

F_k(x) = a_k + b_k(x-x_k) + c_k(x-x_k)² + d_k(x-x_k)³

F = F₁ на интервале [x₀, x₁]
F = F₂ на интервале [x₁, x₂]
...
F = F_N на интервале [x_N-1, x_N]

При этом, на каждом из отрезков [x_k-1; x_k] коэффициенты полинома a_k, b_k, c_k, d_k разные.

Чтобы узнать эти коэффициенты, кроме условия непрерывности функции на полиномы налагают дополнительные условия, а именно непрерывности первой и второй производной функции F(x), а также равенства вторых производных функции на концах отрезка [x₀; x_N], т.е.

F_k-1(x_k-1) = F_k(x_k-1),
F'_k-1(x_k-1) = F'_k(x_k-1),
F''_k-1(x_k-1) = F''_k(x_k-1),
при k=2,3,..N
F''(x₀)=0, F''(x_N)=0.

Найдем выражения для производных функций F_k
F'_k(x) = b_k + 2c_k(x - x_k) + 3d_k(x-x_k)²
F''_k(x) = 2c_k + 6 d_k(x-x_k)

Подставив их в условия непрерывности получим систему:
a₁ - b₁h₁ + c₁h₁² - d₁h₁³ = y₀
a_k = y_k, k=1,2,..N
a_k-1 = a_k - b_kh_k + c_kh_k² - d_kh_k³, k=2,3...N
b_k-1 = b_k - 2c_kh_k + 3d_kh_k², k=2,3...N
c_k-1 = c_k - 3d_kh_k, k=2,3...N
c₁ - 3d₁h₁ = 0
c_N = 0

Здесь введено обозначение h_k = x_k - x_k-1, k=1,2,...N

Введем еще l_k = (y_k - y_k-1)/h_k, k=1,2,...N, а также c₀ = 0 .

Упомянутую выше систему уравнений можно решить с помощью некоторых преобразований, на которых я не буду останавливаться. При этом используется так называемый метод прогонки.

Вводятся прогоночные коэффициенты

δ₁ = - h₂/(2(h₁+h₂))
λ₁ = 3(l₂ - l₁)/(2(h₁+h₂))
δ_k-1 = - h_k/(2h_k-1+2h_k+h_k-1δ_k-2), k=3,4, ... N
λ_k-1 = (3l_k - 3l_k-1 - h_k-1λ_k-2)/(2h_k-1+2h_k+h_k-1δ_k-2)

Далее следует найти коэффициенты c_k по формулам обратной прогонки

c_k-1 = δ_k-1c_k + λ_k-1, k = N, N-1, N-2, ... 2

После нахождения c_k нужно найти b_k и d_k по формулам

b_k = l_k + (2c_kh_k + h_kc_k-1)/3, k=1,2,...N
d_k = (c_k - c_k-1)/(3h_k), k=1,2,...N

Приведем программу, где реализован этот алгоритм:

//Cubic spline interpolation program
//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)  
//where f(xi) = yi
//and f(x) is cubic function on every [x_k-1, x_k] segment
//and f(x), f'(x), f''(x) are continual
//the result is four columns of cubic polinom coefficients

#include <math.h>
#include <stdio.h>
#include <process.h>
float *x, *y, *h, *l, *delta, *lambda, *c, *d, *b;
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));
   N--;
}
void readmatrix(){
   int i=0;
   //read matrixes a and b from input file
   for(i=0; i<N+1; i++){
       fscanf(InFile, "%f", &x[i]);
       fscanf(InFile, "%f", &y[i]);
   }
}

void allocmatrix(){
   //allocate memory for matrixes
   x = new float[N+1];
   y = new float[N+1];
   h = new float[N+1];
   l = new float[N+1];
   delta = new float[N+1];
   lambda = new float[N+1];
   c = new float[N+1];
   d = new float[N+1];
   b = new float[N+1];
}
void freematrix(){
   delete [] x;
   delete [] y;
   delete [] h;
   delete [] l;
   delete [] delta;
   delete [] lambda;
   delete [] c;
   delete [] d;
   delete [] b;
}

void printresult(){
   int k=0;
   printf("\nA[k]\tB[k]\tC[k]\tD[k]\n");
   for(k=1; k<=N; k++){
       printf("%f\t%f\t%f\t%f\n", y[k], b[k], c[k], d[k]);
   }
}
void testresult(){
   float start = x[0];
   float end = x[N];
   float step = (end - start)/20;
   FILE* OutFile = fopen("test.txt", "wt");
   for(float s = start; s<=end; s+= step){
       //find k, where s in [x_k-1; x_k]
       for(int k=1; k<=N; k++){
	   if(s>=x[k-1] && s<=x[k]){
	       break;
	   }
       }
       float F = y[k] + b[k]*(s-x[k]) + c[k]*pow(s-x[k], 2) + d[k]*pow(s-x[k], 3);
       fprintf(OutFile, "%f\t%f\n", s,  F);
   }
   fclose(OutFile);
}
void cls(){
   for(int i=0; i<25; i++) printf("\n");
}
void main(){
   int k=0;
   cls();
   do{
       printf("\nInput filename: ");
       scanf("%s", filename);
       InFile = fopen(filename, "rt");
   }while(InFile==NULL);
   count_num_lines();
   rewind(InFile);
   allocmatrix();
   readmatrix();
   for(k=1; k<=N; k++){
       h[k] = x[k] - x[k-1];
       if(h[k]==0){
	   printf("\nError, x[%d]=x[%d]\n", k, k-1);
	   return;
       }
       l[k] = (y[k] - y[k-1])/h[k];
   }
   delta[1] = - h[2]/(2*(h[1]+h[2]));
   lambda[1] = 1.5*(l[2] - l[1])/(h[1]+h[2]);
   for(k=3; k<=N; k++){
      delta[k-1] = - h[k]/(2*h[k-1] + 2*h[k] + h[k-1]*delta[k-2]);
      lambda[k-1] = (3*l[k] - 3*l[k-1] - h[k-1]*lambda[k-2]) /
		    (2*h[k-1] + 2*h[k] + h[k-1]*delta[k-2]);
   }
   c[0] = 0;
   c[N] = 0;
   for(k=N; k>=2; k--){
      c[k-1] = delta[k-1]*c[k] + lambda[k-1];
   }
   for(k=1; k<=N; k++){
      d[k] = (c[k] - c[k-1])/(3*h[k]);
      b[k] = l[k] + (2*c[k]*h[k] + h[k]*c[k-1])/3;
   }
   printresult();
   testresult();
   freematrix();
}

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

x0	y0
x1	y1
x2	y2
...
xN	yN

Например:

1	5
2	3
3	2.5
4	2
5	0

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

A[k]	B[k]	C[k]	D[k]
3.000000	-1.250000	1.125000	0.375000
2.500000	-0.125000	0.000000	-0.375000
2.000000	-1.250000	-1.125000	-0.375000
0.000000	-2.375000	0.000000	0.375000

Это и есть решение системы, т.е. набор коэффициентов a_k, b_k. c_k, d_k на четырех отрезках.

Кроме того, программа создаст файл test.txt в котором запишет подсчитанные значения интерполирующей функции в 20-точках на отрезке [x₀; x_N]. Вы сможете убедиться, что значения интерполирующей функции плавно меняются от точки к точке, и в точках x_k совпадают со значениями y_k.

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

• void count_num_lines() - подсчитывает количество точек, где задана функция
• void allocmatrix() - выделяет память для массивов b, c, d, x, y, delta, lambda
• void readmatrix() - прочитывает из файла координаты x, y точек
• void testresult() - на основании вычисленных коэффициентов a, b, c, d вычисляет значение интерполирующей функции в 20 точках на интервале [x₀; x_N] с равномерным шагом
• void printresult() - распечатывает столбцы коэффициентов a, b, c, d
• void freematrix() - освобождает память, которая была выделена ранее
• void cls() - стирает экран в начале работы программы
• void main() - основная функция из которой последовательно вызываются все вышеперечисленные функции, и проходит процесс вычисления коэффициентов по формулам прямой и обратной прогонки.

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