# slau yacobi

May 8th, 2021
590
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];
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 << x[i] << "\n";
110.     cout << "\n";
111. }
112. void Yacobi(const int& n, const int& N, const double& c, const double& l, const double& a0, const double& a1, const double& K)
113. {
114.     int i;
115.     //Задаем вектор нач. приближения
116.     complex<double>* x0 = new complex<double>[N]; //по умолчанию иниц. нулями
117.     complex<double>* x = parallel_area(x0, n, N, c, l, a0, a1, K);
118.     while (!mensheeps(x, x0, N))
119.     {
120.         for (i = 0; i < N; i++)
121.             x0[i] = x[i];
122.         x = parallel_area(x0, n, N, c, l, a0, a1, K);
123.     }
124.     show(x, N);
125.     delete[]x0;
126.     delete[]x;
127. }
128. int main()
129. {
130.     setlocale(LC_ALL, "RU");
131.     double c, l, a0, a1, K;
132.     cout << "Введите параметры с, l, a0, a1, к cоответственно: \n";
133.     cin >> c >> l >> a0 >> a1 >> K;
134.     int n, N;
135.     cout << "Введите размерность системы N: \n";
136.     cin >> N;
137.     cout << "Задайте количество потоков: \n";
138.     cin >> n;
139.     Yacobi(n, N, c, l, a0, a1, K);
140.     /*
141.     //Для проверки Inm,Jnm,gamma
142.     ofstream out("sss.txt");
143.     for (int i = 1; i <= 5; i++)
144.     {
145.         cout << "n = " << i << " gamma = " << gamma(i, K) << "\n";
146.         for (int j = 1; j <= 5; 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