Построение
кубического сплайна функции
План:
вывод расчётных формул;
текст программы;
тестирование.
Текст программы.
#include
#include
#include
#include
#include
#include "mat_vec.h" // классы для работы с матрицами и векторами
#include "progonka.h" // решение системы ур-ний (для 3-х диагональных матриц)
#include "funct.h" // второстепеннные функции программы (рисование и т.д.)
// "корень" программы
void spline (float step, int dop, int n, double* &x,double* &y,double* &x1,double* &y1) {
int k = 0;
matrica Sp(n, n-1);
for (int i = 1; i <= (n-1); i++) {
Sp(i,n) = 3*(y[i-1] - 2*y[i] + y[i+1])/pow(step,2);
Sp(i,i) = 4;
if (i < (n-1)) Sp(i,i+1) = 1;
if (i > 1) Sp(i,i-1) = 1;
}
float *tmp;
progonka(Sp, tmp); // решение системы уравнений методом прогонки
// (см. файл "progonka.h")
vector a(n),b(n+1),c(n),d(n); // вычисление коэф-тов многочленов
b(1) = 0;
b(n+1) = 0;
for(int index = 0; index < n-1; index++)
b(index+2) = tmp[index];
delete [] tmp;
for (i = 1; i <= n; i++)
{
d(i) = y[i-1];
a(i) = (b(i+1)- b(i))/(-3*step);
c(i) = (y[i] - d(i) - pow(step,2)*b(i) + pow(step,3)*a(i) )/(-step);
}
i=0;
//построение графика сплайна при помощи полученный коэф-тов (см. выше)
for (i=0; i < n; i++)
for (int j=0; j < dop; j++)
{
x1[k] = x[i] + j*step / (dop);
y1[k] = pow((x[i]-x1[k]),3)*a(i+1)
+ pow((x[i]-x1[k]),2)*b(i+1) + (x[i]-x1[k])*c(i+1)+d(i+1);
k++;
}
x1[n*dop] = x[n];
y1[n*dop] = y[n];
}
void main() {
int n,dop; double step;
cout << "Введите количество интервалов: "; cin >> n;
cout << "Введите количество доп. т. на интервале: "; cin >> dop;
cout << "Введите шаг интервала: "; cin >> step;
dop++;
double *x,*y, *x1,*y1;
initial(x,y,x1,y1,n,dop);
int i = 0; while (i < (n+1)) { // расчёт первоначальных значений функции
x[i] = (i-n/2)*(step);
y[i] = cos(x[i])*pow(x[i],2);
i++;
}
spline (step, dop, n, x,y,x1,y1);
init(); interface(n, dop,x,y,x1,y1);
delete x,y,x1,y1;
closegraph();
}
#ifndef __FUNCT_H
#define __FUNCT_H
#include
// инициализация графики
void init() {
int D,M; D = DETECT; M = 5;
initgraph(&D,&M,"");
}
// рисование графика функции и сплайна
void paint(int Fx,int Fy,int key,int n, int dop, double* &x,double* &y,double* &x1,double* &y1) {
int i = 0, a, b;
a = getmaxx()/2; b = getmaxy()/2;
setfillstyle(0,0); bar(0,0,a*2+1,b*2+1); setcolor(5);
if ((key == 1) || (key == 3))
while ( i < n ) {
line(x[i]*Fx + a, -y[i]*Fy + b, x[i+1]*Fx + a, -y[i+1]*Fy + b);
i = i++;
}
if ((key == 2) || ( key == 3)) {
i = 0;
setcolor(3);
while ( i < n*dop ) {
line(x1[i]*Fx + a, -y1[i]*Fy + b, x1[i+1]*Fx + a, -y1[i+1]*Fy + b);
i = i++;
}
}
setcolor(10); line(getmaxx()/2,0,getmaxx()/2,getmaxy());
line(0,getmaxy()/2,getmaxx(),getmaxy()/2);
}
// функция для приближения (удаления) и масштабирования по осям графиков
void interface(int n, int dop, double* &x, double* &y,double* &x1, double* &y1) {
int c=16, z=16;
char key='0';
while (key != 27) {
if (key == 75) c = c+4;
if (key == 72) z = z+4;
if (key == 77) c = c-4;
if (key == 80) z = z-4;
if (key == 45) { z = z-4; c = c-4; }
if (key == 61) { z = z+4; c = c+4; }
if (c < 0) c = 0;
if (z < 0) z = 0;
if (key == 's') paint(c,z,2,n,dop,x,y,x1,y1);
else if (key == 'f') paint(c,z,1,n,dop,x,y,x1,y1);
else paint(c,z,3,n,dop,x,y,x1,y1);
key = getch();
}
}
// Инициализация динамических массивов
void initial (double* &x,double* &y,double* &x1,double* &y1, int n, int dop) {
x = new double [n+1];
y = new double[n+1];
for (int i = 0 ; i < (n+1);i++) {
y[i] = 0;
x[i] = 0; }
x1 = new double[n*dop+1];
y1 = new double[n*dop+1];
for ( i = 0 ; i < (n*dop+1);i++) {
x1[i] = 0;
y1[i] = 0; }
}
#endif
#ifndef __MAT_VEC_H
#define __MAT_VEC_H
#include
#include
// класс матриц
class matrica {
public:
const int Column, String; //кол-во столбцов и строк матрицы
matrica(int column, int string);
~matrica();
private:
float **WW;
matrica(const matrica& rhs);
matrica& operator=(const matrica& rhs);
public:
float& operator()(int i, int j);
friend ostream& operator<<(ostream& out, const matrica& matr);
friend istream& operator>>(istream& in, const matrica& matr);
};
// конструктор
matrica :: matrica(int column, int string) : Column(column), String(string) {
WW = new float*[string];
if(!WW) {
cout << "\n !!! Не хватает памяти конструктору matrica\n";
exit(EXIT_FAILURE);
}
for(int i = 0; i < string; i++) {
WW[i] = new float[column];
if(!WW[i]) {
cout << "\n !!! Не хватает памяти конструктору matrica\n";
exit(EXIT_FAILURE);
}
for(int j = 0; j < column; j++)
WW[i][j] = 0;
}
}
// деструктор
matrica :: ~matrica() {
for(int i = 0; i < String; i++)
delete [] WW[i];
delete [] WW;
}
// операция доступа к элементу
float& matrica :: operator()(int i, int j) {
if((i > 0) && (i <= String) && (j > 0) && (j <= Column))
return WW[i - 1][j - 1];
else {
cout << "\n Ошибка доступа к элементу (" << i << ", " << j << ") ! \n";
exit(EXIT_FAILURE);
}
}
// вывод матрицы в поток
ostream& operator<<(ostream& out, matrica& WW) {
for(int i = 1; i <= WW.String; i++) {
for(int j = 1; j <= WW.Column; j++)
out << WW(i, j) << " ";
out << endl;
}
return out << "";
}
// ввод матрицы из потока
istream& operator>>(istream& in, matrica& WW) {
for(int i = 1; i <= WW.String; i++)
for(int j = 1; j <= WW.Column; j++)
in >> WW(i, j);
return in;
}
// класс векторов
class vector {
public:
vector(int column);
~vector();
const int Column; // кол-во элементов вектора
private:
float *vect;
vector(const vector& rhs);
vector& operator=(const vector& rhs);
public:
float& operator()(int i);
friend ostream& operator<<(ostream& out, const vector& vec);
friend istream& operator>>(istream& in, const vector& vec);
};
// кнструктор vector
vector :: vector(int column) : Column(column) {
vect = new float[column];
if(!vect) {
cout << endl << "\n !!!Не хватает памяти конструктору vector! \n";
exit(EXIT_FAILURE);
}
for(int i = 0; i < Column; i++)
vect[i] = 0;
}
// деструктор
vector :: ~vector() {
delete [] vect;
}
// операция доступа к эелементу
float& vector :: operator()(int i) {
if((i > 0) && (i <= Column))
return vect[i - 1];
else {
cout << "\n !!!Ошибка доступа к элементу вектора - " << i;
exit(EXIT_FAILURE);
}
}
// вывод вектора в поток
ostream& operator << (ostream& out, vector& vec) {
for(int i = 1; i <= vec.Column; i++)
out << vec(i) << ' ';
return out << endl;
}
// ввод вектора из потока
istream& operator>>(istream& in, vector& vec) {
for(int i = 1; i <= vec.Column; i++)
in >> vec(i);
return in;
}
#endif
#ifndef __PROGONKA_H
#define __PROGONKA_H
#include "mat_vec.h"
int progonka(matrica &mat, float* &x) {
x = new float[mat.String];
if(!x)
return 0;
int i, y = mat.Column, n = mat.String; vector h(n), d(n);
d(1) = - mat(1, 2) / mat(1, 1);
h(1) = mat(1, y) / mat(1, 1);
for(i = 2; i <= n - 1; i++) {
d(i) = mat(i, i+1) / (mat(i, i-1) * d(i-1) - mat(i, i));
h(i) =(mat(i, y)-mat(i,i-1) * h(i-1))/(mat(i, i-1) * d(i-1) + mat(i, i));
}
h(n) =(mat(n, y)-mat(n,n-1) * h(n-1))/(mat(n, n-1) * d(n-1) + mat(n, n));
x[n-1] = h(n); for ( i=n - 1; i >= 1; i--)
x[i - 1] = d(i) * x[i] + h(i);
return 1;
}
#endif
Тестирование:
Зеленым
цветом – график функции
построенный в пределе от –5 до 5, с шагом
= 1.
Красным цветом – график сплайна, полученный при интерполировании исходного графика, причём дополнительно построено всего 3 точки на каждом интервале.