Программная реализация класса аппроксимирующих функций

Класс аппроксимирующих функций должен, по-видимому, содержать в качестве членов-данных два вектора одинаковой размерности для хранения последовательности узлов и значений заданной функции в этих узлах. Набор функций-членов, помимо обязательного набора конструкторов, должен включать методы реализации рассмотренных алгоритмов и, возможно, подпрограммы для интерфейса пользователя. Ниже приведено содержание включаемого файла с определением класса аппроксимаций, содержащего указанные компоненты с сопровождающими комментариями.

#include "polynom.h"

#include "matrix.h"

/*Класс полиномиальных аппроксимаций */

template <class YourOwnFloatType>

class Approximate

{

/*Для решения задачи аппроксимации нам понадобится массивы узлов и значений функции в узлах, т.е. объекты класса vector */

vector<YourOwnFloatType> x, y;/* Узлы и значения в узлах*/

public:/* Общедоступные методы */

/*Основной конструктор в качестве параметра принимает два вектора - вектор узлов и вектор значений в узлах */

Approximate(vector<YourOwnFloatType> _x, vector <YourOwnFloatType> _y): x(_x), y(_y)

{

/*Проверка размерностей на корректность */

if(x.getm()!=y.getm())

throw xmsg("Количество узлов не совпадает с количеством значений в них");

/*Наверное, в этом случае задача аппроксимации теряет смысл: */

if(x.getm()<2)

throw xmsg("Слишком мало точек");

}

/*Конструктор копирования */

Approximate(Approximate &_): x(_.x), y(_.y) {} /*Результатом решения должен быть, естественно, объект класса polynom */

//Прототипы методов класса аппроксимаций

//Интерполяция каноническим полиномом

polynom<YourOwnFloatType> classic(),

//Интерполяция полиномом Ньютона

newton();

/*Интерполяция кубическими сплайнами. Возвращаемое значение - массив полиномов, организуемый динамически. По окончании работы память из под него в обязательном порядке необходимо освободить */

polynom<YourOwnFloatType>* spline();

/*Аппроксимация функции одной переменной по методу наименьших квадратов при степенном базисе. Аргумент - порядок полинома */

polynom<YourOwnFloatType> MNK(int);

};

/*Теперь определения методов класса полиномиальных аппроксимаций */

//Интерполяция каноническим полиномом

template <class YourOwnFloatType>

polynom<YourOwnFloatType> Approximate<YourOwnFloatType>::classic()

{

/*Матрицы для хранения левой и правой частей системы уравнений */

matrix<YourOwnFloatType> mtr(x.getm(),x.getm()), f(x.getm(),1);

//Формирование матрицы СЛАУ

for(int i=0;i<x.getm();i++)

{

/*Правая часть системы - значения функции в узлах */

f[i][0]=y[i];

/*Левая часть системы - вектор степеней узлов */

for(int j=0;j<x.getm();j++)

mtr[i][j]=pow(x[i],j);

}

matrix<YourOwnFloatType> sol=SLAE_Gauss(mtr,f);

/*Вызываем метод Гаусса из матричного класса для решения СЛАУ. /*Результат решения построенной системы - коэффициенты полинома */

polynom<YourOwnFloatType> res=x.getm();

/*формируем полином из коэффициентов */

for(int i=0;i<x.getm();i++)

res[i]=sol[i][0];

return res;/*Возвращаем результат */

}

//Интерполяция полиномом Ньютона

template <class YourOwnFloatType>

polynom<YourOwnFloatType> Approximate<YourOwnFloatType>::newton()

{

//Формирование коэффициентов полинома

matrix<YourOwnFloatType> table(x.getm(),x.getm());

/*Первый столбец таблицы - значения в узлах, во всех последующих - разделённые разности соответствующего порядка */

for(int i=0;i<x.getm();i++)

table[i][0]=y[i];

for(int j=1;j<x.getm();j++)

for(int i=j;i<x.getm();i++)

table[i][j]=

(table[j-1][j-1]-table[i][j-1])/(x[j-1]-x[i]);

/*В диагонали сформированной матрицы находятся коэффициенты полинома в Ньютоновой форме. Для того, чтобы перейти к канонической записи полинома, проделаем следующее: */

/*Объявим и инициализируем два вспомогательных полинома - полиномиальный Х и полиномиальную единицу */

polynom<YourOwnFloatType> Odin,X=2,p;

Odin[0]=1,X[1]=1; /*инициализируем их */

/*сформируем искомый полином как сумму произведений одночленов соответствующих степеней */

for(int i=0;i<x.getm();i++)

{

polynom<YourOwnFloatType> temp=table[i][i]*Odin;

for(int j=0;j<i;j++)

temp*=(X-x[j]*Odin);

p+=temp;

}

return p;/* возвращаем результат */

}

//Интерполяция кубическими сплайнами.

/*Общее количество сплайнов равно количеству межузловых промежутков, т.е. x.getm()-1. В связи с этим результатом выполнения данной функции будет не один полином, а x.getm()-1*/

template <class YourOwnFloatType>

polynom<YourOwnFloatType>* Approximate<YourOwnFloatType>::spline()

{

/*создаём две матрицы для левой и правой частей СЛАУ, результатом решения которой будут производные второго порядка каждого сплайна */

matrix<YourOwnFloatType> mtr(x.getm(),x.getm()), right(x.getm(),1);

/*0,n-1 - в этих точках вторая производная равна нулю*/

mtr[0][0]=1,right[0][0]=0;

mtr[x.getm()-1][x.getm()-1]=1,

right[x.getm()-1][0]=0;

for(long i=1;i<mtr.getm()-1;i++)

{

YourOwnFloatType hi=x[i+1]-x[i],

him1=x[i]-x[i-1];

mtr[i][i-1]=him1;

mtr[i][i]=2*(him1+hi);

mtr[i][i+1]=hi;

right[i][0]=

6*(-(y[i]-y[i-1])/him1+(y[i+1]-y[i])/hi);

}

/*Находим вторые производные и выражаем через них и исходные данные коэффициенты сплайна в форме Тейлора */

matrix<YourOwnFloatType> d2y=SLAE_Orto(mtr,right);

vector<YourOwnFloatType> a=x.getm()-1;

vector<YourOwnFloatType> b=a,c=a,d=a;

for(long i=0;i<a.getm();i++)

{

YourOwnFloatType hi=x[i+1]-x[i];

a[i]=y[i];

b[i]=(y[i+1]-y[i])/hi-hi*(d2y[i+1][0]+2*

d2y[i][0])/6;

c[i]=d2y[i][0]/2;

d[i]=(d2y[i+1][0]-d2y[i][0])/(6*hi);

}

/*Выделяем память под массив полиномов */

polynom<YourOwnFloatType> *arr=

new polynom <YourOwnFloatType>[a.getm()];

/*Составляем два служебных полинома - полиномиальный Х и полиномиальную единицу */

polynom<YourOwnFloatType> Odin,X=2;

Odin[0]=1,X[1]=1;

/*Формируем полиномы, описывающие элементарные сплайны */

for(int i=0;i<a.getm();i++)

arr[i]=a[i]*Odin+b[i]*(X-x[i]*Odin)+c[i]*

((X-x[i]*Odin)^2)+d[i]*((X-x[i]*Odin)^3);

return arr;/*Возвращаем указатель на заполненный массив - не забудьте его освободить! */

}

//Аппроксимация по методу наименьших квадратов

/*Аргумент - степень полинома, которым будем аппроксимировать экспериментально полученные данные */

template <class YourOwnFloatType>

polynom<YourOwnFloatType> Approximate<YourOwnFloatType>::MNK(int rng)

{

matrix<YourOwnFloatType> mtr(x.getm(),rng+1), /*Прямоугольная матрица измерений */

f(x.getm(),1); //Столбец свободных членов

//Формирование матрицы измерений

for(int i=0;i<mtr.getm();i++)

{

f[i][0]=y[i];

for(int j=0;j<mtr.getn();j++)

mtr[i][j]=pow(x[i],j);

}

//Реализуем формулу a=((mtr*mtrT)^-1)*(mtrT*f)

matrix<YourOwnFloatType> sol=

SLAE_Orto((~mtr)*mtr,(~mtr)*f);

polynom<YourOwnFloatType> a(sol.getm()); //Полином - решение

for(long i=0;i<sol.getm();i++)

a[i]=sol[i][0];

return a;/* Возвращаем полином - решение задачи аппроксимации*/

}

/*Определяем типы данных "действительный полином" и "действительный вектор" */

typedef polynom<double> dpolynom;

typedef vector<double> dvector;

#include <conio.h>

/*Тестирующая программа */

void main()

{

/*Узлы и значения в узлах */

double xdata[5]={1,2,3,4,5}, ydata[5]={0.5,1,2,3,3.5};

dvector x(5,xdata),y(5,ydata);

Approximate<double> test(x,y);/*Объявляем тестовый объект аппроксимационного класса */

dpolynom what=test.classic();/*Канонический полином */

cout<<what<<endl<<what(2)<<endl;

what=test.newton(); /*Ньютоновский полином */

cout<<what<<endl<<what(2)<<endl;

what=test.MNK(1); /*МНК, прямая */

cout<<what<<endl<<what(2)<<endl;

what=test.MNK(2); /*МНК, парабола */

cout<<what<<endl<<what(2)<<endl;

what=test.MNK(3); /*МНК, кубическая парабола */

cout<<what<<endl<<what(2)<<endl;

dpolynom *arr=test.spline(); /*получаем указатель на массив сплайнов */

cout<<endl<<endl;

for(long i=0;i<x.getm()-1;i++) /*выводим сплайны с попутным тестом их в узлах */

cout<<arr[i]<<" "<<arr[i](x[i])

<<" "<<arr[i](x[i+1])<<endl;

delete[] arr;/*освобождаем память из под массива сплайнов */

getch();

}