Каталог :: Математика

Реферат: Построение кубического сплайна функции

ГОСУДАРСТВЕННЫЙ КОМИТЕТ РОССИЙСКОЙ ФЕДЕРАЦИИ
ПО ВЫСШЕМУ И СРЕДНЕСПЕЦИАЛЬНОМУ ОБРАЗОВАНИЮ
              КРАСНОЯРСКИЙ ГОСУДАРСТВЕННЫЙ ТЕХНИЧЕСКИЙ УНИВЕРСИТЕТ              
     

Практическое задание

по компьютерной алгебре тема: построение кубического сплайна функции Выполнил: студент 2 курса ФИВТ группы 27-4 Попов А.В. Проверила: Быкова Е.Г. 1998 г. План: 1) вывод расчётных формул; 2) текст программы; 3) тестирование. Текст программы. #include <iostream.h> #include <fstream.h> #include <conio.h> #include <math.h> #include <dos.h> #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 <graphics.h> // инициализация графики 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 <stdlib.h> #include <iostream.h> // класс матриц 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 точки на каждом интервале.