codisinmyvines

SLAU

May 7th, 2021
581
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1. #include <iostream>
  2. #include <omp.h>
  3. #include <cmath>
  4. #include <complex>
  5. #define pi 3.14159265
  6. #define eps 0.000001
  7. using namespace std;
  8.  
  9. complex<double> gamma(double n, const double& K)
  10. {
  11.     complex<double> res;
  12.     if (K * K >= pi * pi * n * n)
  13.     {
  14.         res = sqrt(K * K - pi * pi * n * n);
  15.     }
  16.     else
  17.     {
  18.         res = 1i * sqrt(pi * pi * n * n - K * K);
  19.     }
  20.     return res;
  21. }
  22. complex<double> Jnm(double n, double m, const double& a0, const double& a1)
  23. {
  24.     complex<double> res;
  25.     if (n == m)
  26.     {
  27.         res = a1 - a0 - ((sin((n + m) * pi * a1) - sin((n + m) * pi * a0)) / ((n + m) * pi));
  28.     }
  29.     else
  30.     {
  31.         res = (sin((n - m) * pi * a1) - sin((n - m) * pi * a0)) / ((n - m) * pi) - ((sin((n + m) * pi * a1) - sin((n + m) * pi * a0)) / ((n + m) * pi));
  32.     }
  33.     return res;
  34. }
  35. complex<double> Inm(double n, double m, const double& a0, const double& a1)
  36. {
  37.     complex<double> res;
  38.     if (n == m)
  39.     {
  40.         res = 1.0 - Jnm(n, m, a0, a1);
  41.     }
  42.     else
  43.     {
  44.         res = -Jnm(n, m, a0, a1);
  45.     }
  46.     return res;
  47. }
  48. complex<double> sum(complex<double>* x0, int k, const int& N, const double& c, const double& l, const double& a0, const double& a1, const double& K)
  49. {
  50.     int i, j;
  51.     complex<double>s;
  52.     complex<double>h;
  53.     for (i = 0; i < N; i++)
  54.     {
  55.         h = (1.0 - exp(2.0 * 1i * c * gamma((double)i + 1, K)));
  56.         for (j = 0; j < N; j++)
  57.         {
  58.             s = s + (x0[i] * h * gamma((double)j + 1, K) * Jnm((double)i + 1, (double)j + 1, a0, a1) * Inm((double)j + 1, (double)k + 1, a0, a1)) / (1.0 - exp(2.0 * 1i * gamma((double)j + 1, K) * c));
  59.         }
  60.     }
  61.     return s;
  62. }
  63. complex<double>* parallel_area(complex<double>* x0, const int& n, const int& N, const double& c, const double& l, const double& a0, const double& a1, const double& K)
  64. {  
  65.    
  66.     int i;
  67.     complex<double>* x = new complex<double>[N];
  68.     #pragma omp parallel num_threads(n)
  69.     {
  70.         #pragma omp for
  71.         for (i = 0; i < N; i++)
  72.         {
  73.             x[i] = (sum(x0, i, N, c, l, a0, a1, K) + (gamma(l, K) * Inm(l, (double)i + 1, a0, a1))) / gamma((double)i + 1, K);
  74.         }
  75.     }
  76.     return x;
  77.     //Последовательно
  78.     /*
  79.     //complex<double>s1, s2, s3;
  80.     for (i = 0; i < N; i++)
  81.     {
  82.         //s1 = sum(x0, i, N, c, l, a0, a1, K);
  83.         //s2 = gamma(l, K) * Inm(l, (double)i + 1, a0, a1);
  84.         //s3 = gamma((double)i + 1, K);
  85.         x[i] = (sum(x0, i, N, c, l, a0, a1, K) + gamma(l, K) * Inm(l, (double)i + 1, a0, a1)) / gamma((double)i + 1, K);
  86.         //x[i] = (s1 + s2) / s3;
  87.     }
  88.     return x;
  89.     */
  90. }
  91. bool mensheeps(complex<double>* x, complex<double>* x0, int N)
  92. {
  93.     double max = abs(x[0] - x0[0]);
  94.     double temp;
  95.     for (int i = 1; i < N; i++)
  96.     {
  97.         temp = abs(x[i] - x0[i]);
  98.         if (temp > max)
  99.             max = temp;
  100.     }
  101.     if (max <= eps)
  102.         return true;
  103.     else return false;
  104. }
  105. void show(complex<double>* x, int N)
  106. {
  107.     for (int i = 0; i < N; i++)
  108.         cout <<abs(x[i])<< "\n";
  109.     cout << "\n";
  110. }
  111. void Yacobi(const int& n, const int& N, const double& c, const double& l, const double& a0, const double& a1, const double& K)
  112. {
  113.     int i;
  114.     //Задаем вектор нач. приближения
  115.     complex<double>* x0 = new complex<double>[N]; //по умолчанию иниц. нулями
  116.     for (i = 0; i < N; i++)
  117.         x0[i] = 1.0;
  118.     complex<double>* x = parallel_area(x0, n, N, c, l, a0, a1, K);
  119.     while (!mensheeps(x, x0, N))
  120.     {
  121.         for (i = 0; i < N; i++)
  122.             x0[i] = x[i];
  123.         x = parallel_area(x0, n, N, c, l, a0, a1, K);
  124.     }
  125.     show(x, N);
  126.     delete[]x0;
  127.     delete[]x;
  128. }
  129. int main()
  130. {
  131.     setlocale(LC_ALL, "RU");
  132.     double c, l, a0, a1, K;
  133.     cout << "Введите параметры с, l, a0, a1, к cоответственно: \n";
  134.     cin >> c >> l >> a0 >> a1 >> K;
  135.     int n, N;
  136.     cout << "Введите размерность системы N: \n";
  137.     cin >> N;
  138.     cout << "Задайте количество потоков: \n";
  139.     cin >> n;
  140.     Yacobi(n, N, c, l, a0, a1, K);
  141.     /*
  142.     ofstream out("den.txt");
  143.     for (int i = 1; i <= 10; i++)
  144.     {
  145.         //cout << "n = " << i << " gamma = " << gamma(i, K) << "\n";
  146.         for (int j = 1; j <= 10; j++)
  147.         {
  148.             out << "n = " << i << " m = " << j << "\n"
  149.                 << "Jnm = " << Jnm(i, j, a0, a1) << " Inm = " << Inm(i, j, a0, a1) << " \n";
  150.         }
  151.     }
  152.     out.close();
  153.     */
  154.     return 0;
  155. }
RAW Paste Data