constk

Interpolation

Dec 15th, 2020
86
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1. #include <iostream>
  2. #include <iomanip>
  3. #include <vector>
  4. #include <algorithm>
  5.  
  6.  
  7. // Нужно для использования метода прогонки при интерполяции сплайнами
  8. #include "D:\#Cpp_Projects\Interpolation\SlaeSolving\SlaeSolver.h"
  9. #include "D:\#Cpp_Projects\Interpolation\SlaeSolving\SquareMatrix.cpp"
  10. #include "D:\#Cpp_Projects\Interpolation\SlaeSolving\TridiagonalMatrix.cpp"
  11. #include "D:\#Cpp_Projects\Interpolation\SlaeSolving\ExtendedMatrix.cpp"
  12. #include "D:\#Cpp_Projects\Interpolation\SlaeSolving\SlaeSolver.cpp"
  13. #include "D:\#Cpp_Projects\Interpolation\SlaeSolving\Tester.cpp"
  14.  
  15.  
  16. using namespace std;
  17.  
  18.  
  19. double sqr(double x) { return x * x; };
  20. double cb(double x) { return x * x * x; };
  21.  
  22.  
  23. // Список доступных в программе методов интерполяции
  24. enum class Method { Lagranz, Splines };
  25.  
  26.  
  27.  
  28. // Класс, описывающий функцию y = f(x), заданную в виде таблицы
  29. class Table
  30. {
  31. private:
  32. vector<double> x;
  33. vector<double> y;
  34. size_t n = 0;
  35.  
  36. public:
  37. Table(vector<double> x, vector<double> y) : x(x), y(y)
  38. {
  39. if (x.size() != y.size())
  40. throw invalid_argument("Incorrect function table initialization!");
  41. this->n = x.size();
  42. }
  43.  
  44. vector<double> getX() { return x; }
  45. vector<double> getY() { return y; }
  46. size_t getN() { return n; }
  47.  
  48. bool argsAreOrederedAscending() // Упорядочен ли столбец иксов? (по возрастанию)
  49. {
  50. for (vector<double>::iterator i = x.begin(); i != x.end() - 1; ++i)
  51. if (!(*(i + 1) >= *i))
  52. return false;
  53. return true;
  54. }
  55. };
  56.  
  57.  
  58.  
  59. // Класс, содержащий методы для интерполяции таблично заданных функций различными способами
  60. class Interpolator
  61. {
  62. public:
  63. static double lagranzMethod(Table f, double x) // Методом Лагранжа
  64. {
  65. if (f.getN() == 0)
  66. throw invalid_argument("Table must not be empty!");
  67.  
  68. const int n = f.getN();
  69. vector<double> X = f.getX();
  70. vector<double> Y = f.getY();
  71.  
  72. if(x > *std::max_element(X.begin(), X.end()) || x < *std::min_element(X.begin(), X.end()))
  73. throw invalid_argument("X must be inside interpolation section!");
  74.  
  75. double p = 1.0; // Произведение элементов главной диагонали таблицы разностей
  76. for (int i = 0; i < n; ++i)
  77. p *= (x - X[i]);
  78.  
  79. double sum = 0.0; // Сумма частных (y / k), где y - i-ое значение функции, а k - произведение элементов i-ой строки таблицы разностей
  80. for (int i = 0; i < n; ++i)
  81. {
  82. double k = 1.0; // Произведение элементов i-ой строки таблицы разностей
  83. for (int j = 0; j < n; ++j)
  84. if (i != j)
  85. k *= (X[i] - X[j]);
  86. k *= x - X[i]; // Домножаем k на элемент главной диагонали таблицы разностей
  87.  
  88. sum += Y[i] / k; // Находим сумму частных
  89. }
  90.  
  91. return p * sum; // Ответ - произведение суммы частных и произведения эл-тов главной диагонали таблицы разностей
  92. }
  93.  
  94. static double splinesMethod(Table f, double x) // Через сплайны
  95. {
  96. if (f.getN() == 0)
  97. throw invalid_argument("Table must not be empty!");
  98.  
  99. if(!f.argsAreOrederedAscending())
  100. throw invalid_argument("Arguments row in table must be ordered ascending!");
  101.  
  102.  
  103. // Получаем данные о табличной функции
  104. const int n = f.getN();
  105. const vector<double> X = f.getX();
  106. const vector<double> Y = f.getY();
  107.  
  108.  
  109. // Вычисляем h i-ые шаги между соседними точками в таблице
  110. vector<double> h(n - 1);
  111. for (int i = 0; i != n - 1; ++i)
  112. h[i] = X[i+1] - X[i];
  113.  
  114. vector<double> s(n); // Наклоны сплайна
  115.  
  116. // Матрица и столбец свободных членов для трёхдиагональной слау
  117. TridiagonalMatrix matrix(n);
  118. vector<double> column(n);
  119.  
  120. // Первая строка для слау
  121. matrix.at(0, 0) = -4.0 / h[0];
  122. matrix.at(0, 1) = -2.0 / h[0];
  123. column.at(0) = -6.0 * (Y[1] - Y[0]) / (h[0] * h[0]);
  124.  
  125. // Последняя строка для слау
  126. matrix.at(n - 1, n - 2) = 2.0 / h[n - 2];
  127. matrix.at(n - 1, n - 1) = 4.0 / h[n - 2];
  128. column.at(n - 1) = 6.0 * (Y[n - 1] - Y[n - 2]) / (h[n - 2] * h[n - 2]);
  129.  
  130. // i-ая строка слау, 1 <= i <= n-2
  131. for (int i = 1; i != n - 1; ++i)
  132. {
  133. matrix.at(i, i - 1) = 1.0 / h[i - 1];
  134. matrix.at(i, i) = 2 * (1.0 / h[i - 1] + 1.0 / h[i]);
  135. matrix.at(i, i + 1) = 1.0 / h[i];
  136. column.at(i) = 3 * ((Y[i] - Y[i - 1]) / (h[i - 1] * h[i - 1]) + (Y[i + 1] - Y[i]) / (h[i] * h[i]));
  137. }
  138.  
  139. // Передаём нашу слау в класс, который её решит
  140. SlaeSolver solver(ExtendedMatrix(matrix.getMatrix(), column));
  141.  
  142. // Вычисляем наклоны сплайна методом прогонки
  143. s = solver.sweepMethod();
  144.  
  145. // Вычисляем многочлен Эрмита для заданного промежутка
  146. for (int i = 1; i < n; ++i)
  147. if (x >= X[i - 1] && x <= X[i])
  148. return Y[i - 1] * sqr(x - X[i]) * (2 * (x - X[i - 1]) + h[i - 1]) / cb(h[i - 1]) +
  149. s[i - 1] * sqr(x - X[i]) * (x - X[i - 1]) / sqr(h[i - 1]) +
  150. Y[i] * sqr(x - X[i - 1]) * (2 * (X[i] - x) + h[i - 1]) / cb(h[i - 1]) +
  151. s[i] * sqr(x - X[i - 1]) * (x - X[i]) / sqr(h[i - 1]);
  152.  
  153. throw invalid_argument("X must be inside interpolation section!");
  154. }
  155. };
  156.  
  157.  
  158.  
  159. // Класс, описывающий интерполяционный многочлен, полученный при интерполяции функции тем или иным способом
  160. class Polynomial
  161. {
  162. private:
  163. Table f;
  164. Method method;
  165.  
  166. public:
  167. Polynomial(Table f, Method method) : f(f)
  168. {
  169. if (f.getN() == 0)
  170. throw invalid_argument("Table must not be empty!");
  171.  
  172. // Для метода сплайнов нужно проверить, что иксы упорядочены по возрастанию
  173. if(method == Method::Splines && !f.argsAreOrederedAscending())
  174. throw invalid_argument("Arguments row in table must be ordered ascending!");
  175.  
  176. this->method = method;
  177. }
  178.  
  179. // Считаем значения P(x) в заданной точке
  180. double at(double x)
  181. {
  182. try
  183. {
  184. switch (method)
  185. {
  186. case Method::Lagranz: // P(x) сформирован методом Лагранжа
  187. return Interpolator::lagranzMethod(this->f, x);
  188. case Method::Splines: // P(x) сформирован через сплайны
  189. return Interpolator::splinesMethod(this->f, x);
  190. default:
  191. throw invalid_argument("You must specify method of interpolation correctly!");
  192. }
  193. }
  194. catch (invalid_argument & e)
  195. {
  196. cerr << e.what() << endl;
  197. }
  198. }
  199. };
  200.  
  201.  
  202.  
  203. int main()
  204. {
  205. vector<double> x = { 0.5, 1.0, 7.0, 11.0, 12.0, 15.0, 20.0}; // Столбец иксов для таблицы
  206. vector<double> y = { 23.4, 2.6, -16.0, -186.5, 58.9, 132.0, 300.0}; // Столбец игреков для таблицы
  207.  
  208. Table f(x, y); // Таблично заданная функция
  209.  
  210. // Многочлен, полученный путём интерполяции функции методом Лагранжа
  211. Polynomial p(f, Method::Lagranz);
  212.  
  213. // Вывод значений многочлена
  214. cout << " ================== Lagranz method ================== " << endl;
  215. const int width = 15;
  216. cout << setw(width) << "P (0) = " << setw(width) << p.at(0.0) << endl;
  217. cout << setw(width) << "P (2) = " << setw(width) << p.at(2.0) << endl;
  218. cout << setw(width) << "P (11.5) = " << setw(width) << p.at(11.5) << endl;
  219. cout << setw(width) << "P (16) = " << setw(width) << p.at(16.0) << endl;
  220. cout << setw(width) << "P (35) = " << setw(width) << p.at(35.0) << endl;
  221. cout << endl << endl;
  222.  
  223.  
  224.  
  225. // Многочлен, полученный путём интерполяции функции через сплайны
  226. Polynomial p1(f, Method::Splines);
  227.  
  228. // Вывод значений многочлена
  229. cout << " ================== Splines method ================== " << endl;
  230. cout << setw(width) << "P (0) = " << setw(width) << p1.at(0.0) << endl;
  231. cout << setw(width) << "P (2) = " << setw(width) << p1.at(2.0) << endl;
  232. cout << setw(width) << "P (11.5) = " << setw(width) << p1.at(11.5) << endl;
  233. cout << setw(width) << "P (16) = " << setw(width) << p1.at(16.0) << endl;
  234. cout << setw(width) << "P (35) = " << setw(width) << p1.at(35.0) << endl;
  235.  
  236.  
  237.  
  238. getchar();
  239. return 0;
  240. }
RAW Paste Data