Побудова кубічного сплайну функції

Текст програми.

#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] = (in/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 = новий 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 !!! Бракує пам'яті конструктору matrican";

exit(EXIT_FAILURE);

}

for(int i = 0; i < string; i++) {

WW[i] = new float[column];

if(!WW[i]) {

cout << "n !!! Бракує пам'яті конструктору matrican";

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 точки на кожному інтервалі.

© Реферат плюс