На главную >>>


Интерполяция функций сплайном третьего порядка.

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

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

Пусть на некотором отрезке в точках x0, x1, x2, ... xN нам известны значения некоторой функции f(x), а именно y0, y1, y2, ... yN.

Требуется построить интерполирующую функцию F(x), такую, что она принимает в указаных точках те же значения, т.е. F(x0) = y0, F(x1) = y1, ... F(xN) = yN.

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

Fk(x) = ak + bk(x-xk) + ck(x-xk)2 + dk(x-xk)3

F = F1 на интервале [x0, x1]
F = F2 на интервале [x1, x2]
...
F = FN на интервале [xN-1, xN]

При этом, на каждом из отрезков [xk-1; xk] коэффициенты полинома ak, bk, ck, dk разные.

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

Fk-1(xk-1) = Fk(xk-1),
F'k-1(xk-1) = F'k(xk-1),
F''k-1(xk-1) = F''k(xk-1),
при k=2,3,..N
F''(x0)=0, F''(xN)=0.

Найдем выражения для производных функций Fk
F'k(x) = bk + 2ck(x - xk) + 3dk(x-xk)2
F''k(x) = 2ck + 6 dk(x-xk)

Подставив их в условия непрерывности получим систему:
a1 - b1h1 + c1h12 - d1h13 = y0
ak = yk, k=1,2,..N
ak-1 = ak - bkhk + ckhk2 - dkhk3, k=2,3...N
bk-1 = bk - 2ckhk + 3dkhk2, k=2,3...N
ck-1 = ck - 3dkhk, k=2,3...N
c1 - 3d1h1 = 0
cN = 0

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

Введем еще lk = (yk - yk-1)/hk, k=1,2,...N, а также c0 = 0 .

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

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

δ1 = - h2/(2(h1+h2))
λ1 = 3(l2 - l1)/(2(h1+h2))
δk-1 = - hk/(2hk-1+2hk+hk-1δk-2), k=3,4, ... N
λk-1 = (3lk - 3lk-1 - hk-1λk-2)/(2hk-1+2hk+hk-1δk-2)

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

ck-1 = δk-1ck + λk-1, k = N, N-1, N-2, ... 2

После нахождения ck нужно найти bk и dk по формулам

bk = lk + (2ckhk + hkck-1)/3, k=1,2,...N
dk = (ck - ck-1)/(3hk), 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-е) файл, где напишите значения xk, yk, построчно через пробел, приблизительно так:
x0y0
x1y1
x2y2
...
xNyN

Например:
1 5
2 3
3 2.5
4 2
5 0

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

A[k] B[k] C[k] D[k]
3.000000-1.2500001.1250000.375000
2.500000-0.1250000.000000-0.375000
2.000000-1.250000-1.125000-0.375000
0.000000-2.3750000.0000000.375000

Это и есть решение системы, т.е. набор коэффициентов ak, bk. ck, dk на четырех отрезках.

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

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

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

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


На главную >>>

Hosted by uCoz