Алгоритм реализован в виде программы на языке С, которую вы можете скачать в разделе программы.
Опишем сначала, что такое интерполяция и как эта задача решена в данном случае.
Пусть на некотором отрезке в точках 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) искать полином PN(x), степени N, удовлетворяющий описаным выше условиям, т.е. который проходит через все заданные точки.
Если искать функцию в виде F(x) = PN(x) = a0+ a1x1+ a2x2+ ... aNxN то эта интерполяция называется интерполяцией каноническим полиномом.
Выбор многочлена степени N основан на том факте, что через N+1 точку проходит единственная кривая степени N. Подставив значения функции yi в точках xi в выражение полинома получим систему N+1 линейных уравнений.
| a0 + | a1x0 + | a2x02 + | ... | aNx0N = y0 |
| a0 + | a1x1 + | a2x12 + | ... | aNx1N = y1 |
| a0 + | a1x2 + | a2x22 + | ... | aNx2N = y2 |
| ... | ||||
| a0 + | a1xN + | a2xN2 + | ... | aNxNN = yN |
В данной системе известны все yi, xik, неизвестны все ai. Число неизвестных равно числу уравнений, а именно N+1. Уравнения являются линейными относительно ai, следовательно их можно решить, а решение существует и единственно (при некоторых условиях на детерминант системы уравнений).
Для решения системы воспользуемся методом Гаусса в котором система уравнений приводится к треугольному виду. Я уже описывал этот метод. Описываемая в данном разделе программа основана на уже упомянутой, поэтому я не буду подробно останавливаться на сущности метода Гаусса, лишь приведу текст программы и кратко опишу методику работы с ней, а также ее структуру, а именно, для чего служит такое большое количество подпрограмм в этой программе.
//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-е) файл, где напишите значения xi, yi, построчно через пробел, приблизительно так:
| 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
Это и есть решение системы уравнений, т.е. набор неизвестных коэффициентов ai. Т.е. по нашим данным был построен классический интерполяционный многочлен
F(x) = 10 - 7x + 2.25x2 - 0.25x3.
Так случайно получилось, что он третьей степени, хотя для пяти точек в общем виде должен получаться многочлен четвертой степени.
Теперь для того, чтобы узнать значение функции в некоторой промежуточной точке, например x=2.5, которой не было в таблице, достаточно подставить значение x=2.5 в формулу:
F(2.5) = 10 - 7*2.5 + 2.25*2.52 - 0.25*2.53 = 2.65525
что является правдоподобной оценкой значения функции.
Коротко опишем, для чего служит такое большое количество подпрограмм в данной программе:
- void count_num_lines() - подсчитывает количество точек, где задана функция
- void allocmatrix() - выделяет память для массивов a, x, y
- void readmatrix() - прочитывает из файла координаты точек и вычисляет xik
- void printmatrix() - распечатывает систему уравнений
- void diagonal() - делает так, чтобы на главной диагонали не было нулей, чтобы не пришлось однажды делить на ноль в процессе приведения системы к треугольному виду
- void testsolve() - подставляет решение в систему и распечатывает рядом то, что получается в левой части уравнения и для сравнения печатает рядом свободные члены, те, что в правой части уравнения; получившиеся два столбца должны совпадать, если уравнения решены правильно
- void printresult() - распечатывает получившийся столбец решений
- void freematrix() - освобождает память, которая была выделена ранее
- void cls() - стирает экран в начале работы программы
- void main() - основная функция из которой последовательно вызываются все вышеперечисленные функции, и проходит процесс приведения системы уравнений к треугольному виду и обратная прогонка.
Надеюсь, что вам стало все понятно.
На главную >>>