Advertisement
Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- #include <iostream>
- #include <math.h>
- #include <cmath>
- #include <vector>
- #include <Windows.h>
- #include <stdlib.h>
- #include <fstream>
- #include <string.h>
- #include <time.h>
- #include <iomanip>
- using namespace std;
- const double PI = 3.141592653589793238463;
- //Для хранения сетки
- typedef struct Node
- {
- double x, y;
- } Node;
- //Задание функции
- typedef double(*functiontype)(double x);
- double Myfunc(double x)
- {
- return cos(x);
- //return x * x;
- //return sin(x);
- }
- functiontype Func = &Myfunc;
- double degree(double x, int n) {
- for (int i = 1; i < n; i++)
- x *= x;
- return x;
- }
- //Приблизительное значение 2 производной
- double approximate_derivative_2(functiontype* f, double x) {
- double delta = 1e-5;
- return (((*f)(x) - 2 * (*f)(x + delta) + (*f)(x + 2 * delta)) / (delta*delta));
- }
- //Приблизительное значение 2 производной
- double approximate_derivative_1(functiontype* f, double x) {
- double delta = 1e-5;
- return (((*f)(x + delta) - (*f)(x)) / delta);
- }
- //Возвращает значение сплайна в точке
- double Spline(double x, int i, double *Array_steps, double *x_massiv, Node* Array)
- {
- double h2, g1, g2, x1, x2, y1, y2;
- h2 = Array_steps[i + 1]; //этот массив с индексацией от 1, размер - CountSegments, т.е. кол-в точек - 1
- g1 = x_massiv[i]; g2 = x_massiv[i + 1];
- x1 = Array[i].x; x2 = Array[i + 1].x;
- y1 = Array[i].y; y2 = Array[i + 1].y;
- double a, b, c;
- a = (x*x*x)*((g2 - g1) / (6 * h2)) +
- (x*x)*((g1*x2 - g2 * x1) / (2 * h2)) +
- (x)*((g2*x1*x1 - g1 * x2*x2) / (2 * h2) + (g1*h2 - g2 * h2) / 6 + (y2 - y1) / h2) +
- ((g1*x2*x2*x2 - g2 * x1*x1*x1) / (6 * h2) + (g2*h2*x1 - g1 * h2*x2) / 6 + (y1*x2 - y2 * x1) / h2);
- b = (y1*(x2 - x) + y2 * (x - x1)) / h2 + (g1*degree((x2 - x), 3) - h2 * h2*(x2 - x) + g2 * degree((x - x1), 3) - h2 * h2*(x - x1)) / 6 * h2;
- //if (abs(a - b) > 0.05)
- // cout << " Сплайн считает по разному:" << a << " " << b;
- return (b);
- // return (b);
- }
- //Равномерная сетка
- void ValueUniformTable(functiontype* f, Node* Array, double Initial, double End, int CountSegments, double *Array_steps_uni)
- {
- int i;
- double h;
- Array[0].x = Initial;
- Array[0].y = (*f)(Initial);
- Array_steps_uni[0] = 0;
- h = abs(End - Initial) / CountSegments;
- for (int i = 1; i <= CountSegments; i++)
- {
- Array[i].x = Array[i - 1].x + h;
- Array[i].y = (*f)(Array[i].x);
- Array_steps_uni[i] = h;
- }
- /*cout << endl << "Равномерная " << endl;
- for (int i = 0; i < CountSegments+1; i++)
- cout << "(" << Array[i].x << ":" << Array[i].y << ")" << endl;
- */
- }
- //Чебышёвская сетка
- void ValueChebyshevTable(functiontype* f, Node* Array, double Initial, double End, int CountSegments, double *Array_steps_cheb)
- { // Создание таблицы Чебышевских значений
- Array_steps_cheb[0] = 0;
- Array[0].x = Initial;
- Array[0].y = (*f)(Array[0].x);
- for (int i = 1; i < CountSegments; i++)
- {
- Array[i].x = -((End + Initial) / 2)
- - ((End - Initial) / 2) * cos(((2 * i + 1) * PI) / (2 * (CountSegments+1)));
- Array[i].y = (*f)(Array[i].x);
- Array_steps_cheb[i] = Array[i].x - Array[i - 1].x;
- }
- Array[CountSegments].x = End;
- Array[CountSegments].y = (*f)(Array[CountSegments].x);
- Array_steps_cheb[CountSegments] = Array[CountSegments].x - Array[CountSegments-1].x;
- /*cout << endl << "Cheb " << endl;
- for (int i = 0; i < CountSegments+1; i++)
- cout << "(" << Array[i].x << ":" << Array[i].y << ")" << endl;*/
- }
- //Неравномерная сетка
- void ValueIrregularTable(functiontype* f, Node* Array, double Initial, double End, int CountSegments, double *Array_steps)
- {
- int i;
- double alpha, h;
- alpha = 2 * PI / CountSegments;
- Array_steps[0] = 0;
- Array[0].x = Initial;
- Array[0].y = (*f)(Array[0].x);
- Array[CountSegments].x = End;
- Array[CountSegments].y = (*f)(Array[CountSegments].x);
- h = (Array[CountSegments].x - Array[0].x) / CountSegments;
- for (i = 1; i < CountSegments; i++) {
- Array_steps[i] = h + (9 * h / 10)*cos(alpha*i);
- Array[i].x = Array[i - 1].x + Array_steps[i];
- Array[i].y = (*f)(Array[i].x);
- }
- Array_steps[CountSegments] = Array[CountSegments].x - Array[CountSegments - 1].x;
- //вывод в консоль
- /*cout << endl << "Неравномерная" << endl;
- for (i = 0; i <= CountSegments; i++)
- cout << "(" << Array[i].x << ":" << Array[i].y << ")" << endl; */
- }
- /* Принимает на вход коэффициенты СЛАУ в виде массива matrix_coeffs - массив вида [[a1,b1,c1,d1],[a2,b2,c2,d2],...,[an,bn,cn,dn].
- Заполняет Двумерный массив прогоночных коэффициентов coeffs_massiv (0 - Кси, 1 - Эта).
- Находит решение в виде массива x_massiv
- */
- void tridiagonal_matrix_algorithm(double** matrix_coeffs, double* Array_steps, Node* Array, int CountSegments, double* x_massiv, double B) {
- int i;
- double denominator;
- double K, E, K_prev = 0, E_prev = 0; //Прогоночные коэффициенты
- double y1, y2, h2;
- int matrix_size = CountSegments - 1;
- double** coeffs_massiv;
- coeffs_massiv = new double*[matrix_size + 1];
- for (i = 0; i <= matrix_size; i++)
- coeffs_massiv[i] = new double[2];
- coeffs_massiv[1][0] = -(matrix_coeffs[0][2])/(matrix_coeffs[0][1]); //вместо мнимых определил К2,Е2 как следствие из них
- coeffs_massiv[1][1] = (matrix_coeffs[0][3])/(matrix_coeffs[0][1]);
- //Прямой ход
- for (i = 2; i <= matrix_size; i++) { //Ищем K_i+1 и E_i+1 на iом шаге //со 2 начал
- denominator = (matrix_coeffs[i - 1][0] * K_prev + matrix_coeffs[i - 1][1]);
- K = -(matrix_coeffs[i - 1][2]) / denominator;
- E = (matrix_coeffs[i - 1][3] - matrix_coeffs[i - 1][0] * E_prev) / denominator;
- K_prev = K;
- E_prev = E;
- coeffs_massiv[i][0] = K; coeffs_massiv[i][1] = E; //лежат со сдвигом из за индексации
- }
- //Обратный ход
- x_massiv[matrix_size] = coeffs_massiv[matrix_size][1]; //Вместо зануления xn+1 определяю xn = e_n+1
- for (i = matrix_size-1; i >= 1; i--) {
- x_massiv[i] = coeffs_massiv[i][0] * x_massiv[i + 1] + coeffs_massiv[i][1];
- }
- //Определяю gn
- h2 = Array_steps[CountSegments]; //вроде так, h0 условно существует и h1 лежит на 1 инд, h_n-1 na matrix size
- y1 = Array[CountSegments - 1].y;
- y2 = Array[CountSegments].y;
- x_massiv[CountSegments] = 6 * (y2 - y1 - B * h2) / (h2*h2) - 2 * x_massiv[CountSegments - 1]; //gn
- /*cout << endl << "Размер матрицы = " << matrix_size << ", Гамм " << matrix_size + 2 << endl;
- for (i = 0; i < matrix_size + 2; i++)
- cout << "g_ " << i << ": " << x_massiv[i] << endl;
- cout << endl << endl;
- for (i = 0; i < matrix_size+1; i++)
- cout << "e,k_i " << i << ": " << coeffs_massiv[i][0] << " " << coeffs_massiv[i][1] << endl;*/
- }
- //Составляем СЛАУ с учётом нач данных и закидываем это всё в один массив
- void get_matrix_coeffs(double Initial, double End, int CountSegments, functiontype* f, double** matrix_coeffs, Node* Array, double *Array_steps, double B) {
- int i;
- double h1, h2, h_prev, d2, g2;
- double x1, x2, y1, y2, y3, denominator;
- int matrix_size = CountSegments - 1;
- for (i = 1; i <= matrix_size; i++) { //в системе i=1,..n-1, g0,gn - известны
- h1 = Array_steps[i]; //i
- h2 = Array_steps[i + 1];; //i+1
- h_prev = h2;
- matrix_coeffs[i - 1][0] = h1;
- matrix_coeffs[i - 1][1] = 2 * (h1 + h2);
- matrix_coeffs[i - 1][2] = h2;
- matrix_coeffs[i - 1][3] = 6 * ((Array[i + 1].y - Array[i].y) / h2 - (Array[i].y - Array[i - 1].y) / h1);
- }
- matrix_coeffs[0][0] = 0; //а1=0
- matrix_coeffs[matrix_size - 1][2] = 0; //cn=0
- /*cout << "(" << matrix_coeffs[0][1] << ")" << " " << "(" << matrix_coeffs[0][2] << ")" << " " << "(" << "0" << ")" << endl;
- for (i = 1; i < matrix_size - 1; i++) {
- cout << "(" << matrix_coeffs[i][0] << ")" << " " << "(" << matrix_coeffs[i][1] << ")" << " " << "(" << matrix_coeffs[i][2] << ")" << " " << endl;
- }
- cout << "(" << "0" << ")" << " " << "(" << matrix_coeffs[matrix_size - 1][0] << ")" << " " << "(" << matrix_coeffs[matrix_size - 1][1] << ")" << " " << endl; */
- }
- /*
- Строит график функции на всём интервале с заданным количеством точек Count_dots
- */
- void orig_table_in_file(functiontype* f, int Count_dots, double Initial, double End) {
- int k, i;
- double x_value, step;
- step = (End - Initial) / (Count_dots - 1);
- x_value = Initial - step;
- ofstream fout("D:/original_graphic.txt");
- for (i = 0; i < Count_dots; i++) {
- x_value += step;
- fout << x_value << " ";
- fout << (*f)(x_value) << endl;
- }
- fout.close();
- }
- //Используя функцию Spline записываем в файл таблицу значений спалйна на всём интервале в Count_dots точках
- void spline_table_in_file(Node* Array, int Count_dots, int Count_Segments, double Initial, double End, double *Array_steps, double *x_massiv) { // Функция берет таблицу иксов и игриков, количество точек в которых считаем знач полинома (мб убрать) и коэфф полинома,
- int i = 0;
- double step = (End - Initial) / (Count_dots - 1), y_value, x_value;
- ofstream fout("D:/spline_graphic.txt");
- x_value = Initial;
- for (i = 0; i < Count_Segments; i++) {
- while (x_value < Array[i + 1].x) {
- fout << x_value << " ";
- y_value = Spline(x_value, i, Array_steps, x_massiv, Array);
- fout << y_value << endl;
- x_value += step;
- }
- }
- fout << End << " ";
- fout << Spline(End, Count_Segments - 1, Array_steps, x_massiv, Array) << endl;
- fout.close();
- }
- int main()
- {
- setlocale(LC_ALL, "RUS");
- functiontype Func = &Myfunc;
- int CountSegments, i, Countdots = 5000;
- /*---------------------------------Входные данные-------------------------------------*/
- double Initial = -5, End = 5, B, A;
- cout << "Введите N: ";
- //cin >> CountSegments;
- CountSegments = 15000;
- cout << "Точек: " << CountSegments + 1 << endl << endl;
- cout << "Введите A,B: " << endl << endl;
- // cin >> A >> B;
- A = approximate_derivative_2(&Func, Initial); B = approximate_derivative_1(&Func, End);
- /*---------------------------------Сетки-------------------------------------*/
- double *Array_steps_irreg = new double[CountSegments + 1];
- double *Array_steps_cheb = new double[CountSegments + 1];
- double *Array_steps_uni = new double[CountSegments + 1];
- //Равномерная сетка
- Node* ArrayUniformNodes = new Node[CountSegments + 1];
- ValueUniformTable(&Func, ArrayUniformNodes, Initial, End, CountSegments, Array_steps_uni);
- //Чебышёв
- Node* ArrayChebyshevNodes = new Node[CountSegments + 1];
- ValueChebyshevTable(&Func, ArrayChebyshevNodes, Initial, End, CountSegments, Array_steps_cheb);
- //Неравномерная сетка
- Node* ArrayIrregularNodes = new Node[CountSegments + 1];
- ValueIrregularTable(&Func, ArrayIrregularNodes, Initial, End, CountSegments, Array_steps_irreg);
- //Шаги
- /*cout << endl;
- for (i = 0; i <= CountSegments; i++) {
- cout << "h_" << i << ": " << Array_steps_cheb[i] << endl;
- }*/
- /*---------------------------------СЛАУ - matrix_coeffs -------------------------------------*/
- double** matrix_coeffs;
- matrix_coeffs = new double*[CountSegments - 1];
- for (i = 0; i < CountSegments - 1; i++)
- matrix_coeffs[i] = new double[3];
- cout << endl << "Полученная матрица: " << endl;
- get_matrix_coeffs(Initial, End, CountSegments, &Func, matrix_coeffs, ArrayChebyshevNodes, Array_steps_cheb, B);// ДЛЯ ЧЕБЫШЕВА ПРОВЕРИТЬ ВСЕ МАССИВЫ ЧТОБЫ НОРМ ПО НЕМУ ЗАДАВАЛИСЬ (ЭРЭЙ СТЕПС)
- /*cout << "Правая часть: " << endl << endl;
- for (i = 0; i < CountSegments - 1; i++) {
- cout << "d_" << i << ": " << matrix_coeffs[i][3] << endl;
- }*/
- /*---------------------------------Прогонка - запись в x_massiv -------------------------------------*/
- //Получение ответа в x_massiv
- double* x_massiv = new double[CountSegments + 1];
- x_massiv[0] = A; //g0
- tridiagonal_matrix_algorithm(matrix_coeffs, Array_steps_cheb, ArrayChebyshevNodes, CountSegments, x_massiv, B);
- /*---------------------------------Графики -------------------------------------*/
- //Вывод интерполируемой функции
- orig_table_in_file(&Func, Countdots, Initial, End);
- //Вывод сплайна на экран
- spline_table_in_file(ArrayChebyshevNodes, Countdots, CountSegments, Initial, End, Array_steps_cheb, x_massiv);
- cout << endl;
- system("pause");
- return 0;
- }
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement