Advertisement
vatman

for_gui_mmn

May 27th, 2024
634
0
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
C++ 13.03 KB | None | 0 0
  1. #include <cmath>
  2. #include <fstream>
  3. #include <iomanip>
  4. #include <iostream>
  5. #include <vector>
  6.  
  7. extern "C" {
  8. double f_test(double x, double y) {
  9.   double T = M_PI * x * y;
  10.   double sin2T = std::pow(std::sin(T), 2);
  11.   double cos2T = std::pow(std::cos(T), 2);
  12.   double y2_x2 = y * y + x * x;
  13.   return -2 * M_PI * std::exp(std::pow(std::sin(T), 2)) *
  14.          (2 * M_PI * y2_x2 * cos2T * sin2T - M_PI * y2_x2 * sin2T +
  15.           M_PI * y2_x2 * cos2T);
  16. }
  17. }
  18.  
  19. extern "C" {
  20. double f_main(double x, double y) {
  21.   return -std::pow(std::sin(M_PI * y * x), 2);
  22. }
  23. }
  24. extern "C" {
  25. double u(double x, double y) {
  26.   return std::exp(std::pow(std::sin(M_PI * x * y), 2));
  27. }
  28. }
  29. extern "C" {
  30. void Write_test(std::vector<std::vector<double>> v, const int n, const int m,
  31.                 const int S, const double eps_max, const double r_max,
  32.                 const double z_max) // функция записи массива в бинарный файл
  33. {
  34.   std::ofstream in("/home/syatov430/VAZHNO/NM_labs_3kurs_2semestr/NM_lab2/lab2/"
  35.                    "data_test.txt"); // подключение текстового файла для записи
  36.                                      // в бинарном режиме
  37.   for (int i = 0; i < n; i++) {
  38.     for (int j = 0; j < m; j++) {
  39.       if (j != m - 1)
  40.         in << v[i][j] << ", "; // запись в файл поэлементно
  41.       else
  42.         in << v[i][j] << "\n";
  43.     }
  44.   }
  45.   in.close(); // закрываем файл
  46.  
  47.   std::ofstream in1("/home/syatov430/VAZHNO/NM_labs_3kurs_2semestr/NM_lab2/"
  48.                     "lab2/z_it.txt"); // подключение текстового файла для записи
  49.                                       // в бинарном режиме
  50.   for (int i = 0; i < n; i++) {
  51.     for (int j = 0; j < m; j++) {
  52.       if (j != m - 1)
  53.         in1 << fabs(v[i][j] - u(i * (1 / (n - 1)), j * (1 / (m - 1))))
  54.             << ", "; // запись в файл поэлементно
  55.       else
  56.         in1 << fabs(v[i][j] - u(i * (1 / (n - 1)), j * (1 / (m - 1)))) << "\n";
  57.     }
  58.   }
  59.   in1.close(); // закрываем файл
  60.  
  61.   std::ofstream in2("/home/syatov430/VAZHNO/NM_labs_3kurs_2semestr/NM_lab2/"
  62.                     "lab2/param_test.txt"); // подключение текстового файла для
  63.                                             // записи в бинарном режиме
  64.   in2 << S << ", " << eps_max << ", " << r_max << ", " << z_max;
  65.   in2.close(); // закрываем файл
  66. }
  67. }
  68.  
  69. extern "C" {
  70. void Write_main_v(std::vector<std::vector<double>> v, const int n, const int m,
  71.                   const double v_v2, const int n2,
  72.                   const int m2) // функция записи массива в бинарный файл
  73. {
  74.   std::ofstream in("/home/syatov430/VAZHNO/NM_labs_3kurs_2semestr/NM_lab2/lab2/"
  75.                    "data_main.txt"); // подключение текстового файла для записи
  76.                                      // в бинарном режиме
  77.   for (int i = 0; i < n; i++) {
  78.     for (int j = 0; j < m; j++) {
  79.       if (j != m - 1)
  80.         in << v[i][j] << ", "; // запись в файл поэлементно
  81.       else
  82.         in << v[i][j] << "\n";
  83.     }
  84.   }
  85.   in.close(); // закрываем файл
  86.   std::ofstream in2(
  87.       "/home/syatov430/VAZHNO/NM_labs_3kurs_2semestr/NM_lab2/lab2/"
  88.       "param_main_v2.txt"); // подключение текстового файла для записи в
  89.                             // бинарном режиме
  90.   in2 << v_v2 << n2 << m2;
  91.   in2.close(); // закрываем файл
  92. }
  93. }
  94.  
  95. void Write_main_param(
  96.     const int n, const int m, const int S, const double eps_max,
  97.     const double r_max) // функция записи массива в бинарный файл
  98. {
  99.  
  100.   std::ofstream in2("/home/syatov430/VAZHNO/NM_labs_3kurs_2semestr/NM_lab2/"
  101.                     "lab2/param_main.txt"); // подключение текстового файла для
  102.                                             // записи в бинарном режиме
  103.   in2 << n << m << S << ", " << eps_max << ", " << r_max;
  104.   in2.close(); // закрываем файл
  105. }
  106. extern "C" {
  107. int MinimalResiduals_test(const int nmax = 1000, const double _eps = 0.0000001,
  108.                           const int _n = 3, const int _m = 3) {
  109.   int Nmax = nmax; // максимальное число итераций (не менее 1)
  110.   std::vector<std::vector<double>> v; // сеточная функция ѵ (x, y)
  111.   v.resize(_n + 1);
  112.   for (size_t i = 0; i < _n + 1; i++) {
  113.     v[i].resize(_m + 1);
  114.     for (size_t j = 0; j < _m + 1; j++) {
  115.       v[i][j] = 0;
  116.     }
  117.   }
  118.   std::vector<std::vector<double>> r_vec; // вектор невязки
  119.   r_vec.resize(_n + 1);
  120.   for (size_t i = 0; i < _n + 1; i++) {
  121.     r_vec[i].resize(_m + 1);
  122.     for (size_t j = 0; j < _m + 1; j++) {
  123.       r_vec[i][j] = 0;
  124.     }
  125.   }
  126.   std::vector<std::vector<double>> ar; // вектор невязки
  127.   ar.resize(_n + 1);
  128.   for (size_t i = 0; i < _n + 1; i++) {
  129.     ar[i].resize(_m + 1);
  130.     for (size_t j = 0; j < _m + 1; j++) {
  131.       ar[i][j] = 0;
  132.     }
  133.   }
  134.   double r1 = 0;
  135.   int S = 0;         // счетчик итераций
  136.   double eps = _eps; // заданная точность
  137.   double eps_max = 0; // точность, достигнутая на текущей итерации
  138.   double eps_cur = 0; // для подсчета точности на текущей итерации
  139.   double a2, k2, h2; // ненулевые элементы матрицы (-А)
  140.   const int n = _n, m = _m; // размерность сетки
  141.   double r_max = 0;
  142.   double a = 0;
  143.   double b = 1;
  144.   double c = 0;
  145.   double d = 1; // границы области определения уравнения
  146.   int i, j; // индексы
  147.   double r = 0;
  148.   double v_old; // старое значение преобразуемой компоненты вектора
  149.   double v_new; // новое значение преобразуемой компоненты вектора и
  150.   bool flag = false; // условие остановки
  151.   double h = ((b - a) / n);
  152.   double k = ((d - c) / m);
  153.   double z = 0;
  154.   double z_max = 0;
  155.   double tau_s = 0;
  156.   double num = 0; // tau=num/denominator
  157.   double denominator = 1;
  158.   h2 = -std::pow((n / (b - a)), 2);
  159.   k2 = -std::pow((m / (d - c)), 2);
  160.   a2 = -2 * (h2 + k2);
  161.   for (int i = 0; i < n + 1; i++) {
  162.     v[i][0] = u(a + i * h, a);
  163.     v[i][m] = u(a + i * h, b);
  164.   }
  165.   for (int j = 0; j < m + 1; j++) {
  166.     v[0][j] = u(c, c + j * k);
  167.     v[n][j] = u(d, c + j * k);
  168.   }
  169.   do {
  170.     z = 0;
  171.     z_max = 0;
  172.     eps_max = 0;
  173.     r_max = -20;
  174.     num = 0;
  175.     denominator = 0;
  176.     for (j = 1; j < m; j++) {
  177.       for (i = n - 1; i > 0; i--) {
  178.         r_vec[i][j] = f_test(a + i * h, c + j * k) -
  179.                       (h2 * (v[i + 1][j] + v[i - 1][j]) +
  180.                        k2 * (v[i][j + 1] + v[i][j - 1]) + v[i][j] * a2);
  181.       }
  182.     }
  183.     for (j = 1; j < m; j++) {
  184.       for (i = n - 1; i > 0; i--) {
  185.         ar[i][j] =
  186.             (h2 * (r_vec[i + 1][j] + r_vec[i - 1][j]) +
  187.              k2 * (r_vec[i][j + 1] + r_vec[i][j - 1]) + r_vec[i][j] * a2);
  188.         num += ar[i][j] * r_vec[i][j];
  189.         denominator += ar[i][j] * ar[i][j];
  190.       }
  191.     }
  192.     tau_s = num / denominator;
  193.     for (j = 1; j < m; j++) {
  194.       for (i = n - 1; i > 0; i--) {
  195.         v_old = v[i][j];
  196.         v_new = v_old + tau_s * r_vec[i][j];
  197.         eps_cur = std::fabs(v_old - v_new);
  198.         if (eps_cur > eps_max) {
  199.           eps_max = eps_cur;
  200.         }
  201.         v[i][j] = v_new;
  202.         r = fabs(f_test(a + i * h, c + j * k) -
  203.                  (h2 * (v[i + 1][j] + v[i - 1][j]) +
  204.                   k2 * (v[i][j + 1] + v[i][j - 1]) + v_old * a2));
  205.         if (r > r_max)
  206.           r_max = r;
  207.         z = u(a + i * h, c + j * k) - v[i][j];
  208.         if (z > z_max)
  209.           z_max = z;
  210.       }
  211.     }
  212.  
  213.     S = S + 1;
  214.     if ((eps_max <= eps) or (S >= Nmax)) {
  215.       flag = true;
  216.     }
  217.   } while (!flag);
  218.   Write_test(v, n + 1, m + 1, S, eps_max, r_max, z_max);
  219.   return 0;
  220. }
  221. }
  222.  
  223. extern "C" {
  224. std::vector<std::vector<double>>
  225. MinimalResiduals_main(const int nmax = 1000, const double _eps = 0.0000001,
  226.                       const int _n = 3, const int _m = 3) {
  227.   int Nmax = nmax; // максимальное число итераций (не менее 1)
  228.   std::vector<std::vector<double>> v; // сеточная функция ѵ (x, y)
  229.   v.resize(_n + 1);
  230.   for (size_t i = 0; i < _n + 1; i++) {
  231.     v[i].resize(_m + 1);
  232.     for (size_t j = 0; j < _m + 1; j++) {
  233.       v[i][j] = 0;
  234.     }
  235.   }
  236.   std::vector<std::vector<double>> v2; // сеточная функция ѵ (x, y)
  237.   v2.resize(2 * _n + 1);
  238.   for (size_t i = 0; i < 2 * _n + 1; i++) {
  239.     v2[i].resize(2 * _m + 1);
  240.     for (size_t j = 0; j < 2 * _m + 1; j++) {
  241.       v2[i][j] = 0;
  242.     }
  243.   }
  244.   std::vector<std::vector<double>> r_vec; // вектор невязки
  245.   r_vec.resize(_n + 1);
  246.   for (size_t i = 0; i < _n + 1; i++) {
  247.     r_vec[i].resize(_m + 1);
  248.     for (size_t j = 0; j < _m + 1; j++) {
  249.       r_vec[i][j] = 0;
  250.     }
  251.   }
  252.   std::vector<std::vector<double>> ar; // вектор невязки
  253.   ar.resize(_n + 1);
  254.   for (size_t i = 0; i < _n + 1; i++) {
  255.     ar[i].resize(_m + 1);
  256.     for (size_t j = 0; j < _m + 1; j++) {
  257.       ar[i][j] = 0;
  258.     }
  259.   }
  260.   double r1 = 0;
  261.   int S = 0;         // счетчик итераций
  262.   double eps = _eps; // заданная точность
  263.   double eps_max = 0; // точность, достигнутая на текущей итерации
  264.   double eps_cur = 0; // для подсчета точности на текущей итерации
  265.   double a2, k2, h2; // ненулевые элементы матрицы (-А)
  266.   const int n = _n, m = _m; // размерность сетки
  267.   double r_max = 0;
  268.   double a = 0;
  269.   double b = 1;
  270.   double c = 0;
  271.   double d = 1; // границы области определения уравнения
  272.   int i, j; // индексы
  273.   double r = 0;
  274.   double v_old; // старое значение преобразуемой компоненты вектора
  275.   double v_new; // новое значение преобразуемой компоненты вектора и
  276.   bool flag = false; // условие остановки
  277.   double h = ((b - a) / n);
  278.   double k = ((d - c) / m);
  279.   double z = 0;
  280.   double z_max = 0;
  281.   double tau_s = 0;
  282.   double num = 0; // tau=num/denominator
  283.   double denominator = 1;
  284.   h2 = -std::pow((n / (b - a)), 2);
  285.   k2 = -std::pow((m / (d - c)), 2);
  286.   a2 = -2 * (h2 + k2);
  287.   for (int i = 0; i < n + 1; i++) {
  288.     double x = a + i * h;
  289.     v[i][0] = x - x * x;
  290.     v[i][m] = x - x * x;
  291.   }
  292.   for (int j = 0; j < m + 1; j++) {
  293.     double y = c + j * k;
  294.     v[0][j] = std::sin(M_PI * y);
  295.     v[n][j] = std::sin(M_PI * y);
  296.   }
  297.   do {
  298.     z = 0;
  299.     z_max = 0;
  300.     eps_max = 0;
  301.     r_max = -20;
  302.     num = 0;
  303.     denominator = 0;
  304.     for (j = 1; j < m; j++) {
  305.       for (i = n - 1; i > 0; i--) {
  306.         r_vec[i][j] = f_main(a + i * h, c + j * k) -
  307.                       (h2 * (v[i + 1][j] + v[i - 1][j]) +
  308.                        k2 * (v[i][j + 1] + v[i][j - 1]) + v[i][j] * a2);
  309.       }
  310.     }
  311.     for (j = 1; j < m; j++) {
  312.       for (i = n - 1; i > 0; i--) {
  313.         ar[i][j] =
  314.             (h2 * (r_vec[i + 1][j] + r_vec[i - 1][j]) +
  315.              k2 * (r_vec[i][j + 1] + r_vec[i][j - 1]) + r_vec[i][j] * a2);
  316.         num += ar[i][j] * r_vec[i][j];
  317.         denominator += ar[i][j] * ar[i][j];
  318.       }
  319.     }
  320.     tau_s = num / denominator;
  321.     for (j = 1; j < m; j++) {
  322.       for (i = n - 1; i > 0; i--) {
  323.         v_old = v[i][j];
  324.         v_new = v_old + tau_s * r_vec[i][j];
  325.         eps_cur = std::fabs(v_old - v_new);
  326.         if (eps_cur > eps_max) {
  327.           eps_max = eps_cur;
  328.         }
  329.         v[i][j] = v_new;
  330.         r = fabs(f_main(a + i * h, c + j * k) -
  331.                  (h2 * (v[i + 1][j] + v[i - 1][j]) +
  332.                   k2 * (v[i][j + 1] + v[i][j - 1]) + v_old * a2));
  333.         if (r > r_max)
  334.           r_max = r;
  335.       }
  336.     }
  337.  
  338.     S = S + 1;
  339.     if ((eps_max <= eps) or (S >= Nmax)) {
  340.       flag = true;
  341.     }
  342.   } while (!flag);
  343.   Write_main_param(n + 1, m + 1, S, eps_max, r_max);
  344.   return v;
  345. }
  346. }
  347.  
  348. extern "C" {
  349. int do_all() {
  350.  
  351.   int n = 10;
  352.   int m = 10;
  353.   std::vector<std::vector<double>> v_1(n, std::vector<double>(m, 0));
  354.   std::vector<std::vector<double>> v_2(n, std::vector<double>(2 * m, 0));
  355.   v_1 = MinimalResiduals_main(100000, 0.000001, n, m);
  356.   v_2 = MinimalResiduals_main(100000, 0.000001, 2 * n, 2 * m);
  357.   n++;
  358.   m++;
  359.   double maxv_v2 = 0.0;
  360.   double v_v2 = 0.0;
  361.   int n2 = 0;
  362.   int m2 = 0;
  363.   for (int i = n - 1; i >= 0; i--) {
  364.     for (int j = 0; j < m; j++) {
  365.       v_v2 = std::fabs(v_1[i][j] - v_2[i * 2][j * 2]);
  366.       if (v_v2 > maxv_v2) {
  367.         maxv_v2 = v_v2;
  368.         n2 = i;
  369.         m2 = j;
  370.       }
  371.     }
  372.   }
  373.   Write_main_v(v_1, n, m, maxv_v2, n2, m2);
  374.   return 0;
  375. }
  376. }
  377. int main() {
  378.   do_all();
  379.   return 0;
  380. }
  381.  
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement