API tools faq
paste
Advertisement
frustration

squares

frustration
Feb 15th, 2020
217
0
Never
Add comment
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
text 19.42 KB | None | 0 0
raw download clone embed print report
  1. #define _CRT_SECURE_NO_WARNINGS
  2.  
  3. #include <conio.h>
  4. #include <iostream>
  5. #include <cmath>
  6. #include <vector>
  7. #include <fstream>
  8. #include <stdio.h>
  9. #include <string>
  10.  
  11. using namespace std;
  12. const double pi = 3.14159265359;
  13.  
  14. void transposition_polar(vector < vector < double > > &lambda_bb, vector < vector < double > > &lambda_bn,
  15. vector < vector < double > > &lambda_bf, vector < vector < double > >&lambda_fb, vector < vector < double > > &lambda_fn,
  16. vector < vector < double > > &lambda_ff, vector < vector < double > > &lambda_nb, vector < vector < double > > &lambda_nn,
  17. vector < vector < double > > &lambda_nf, vector < double > square_front, vector < double > square_side, vector < double > square_up,
  18. vector <vector < double >> volume,int r, int z) {
  19.  
  20. for (int i = 1; i < r + 1; ++i) {
  21. for (int j = 1; j < z + 1; ++j) {
  22. lambda_bb[i][j] = 0.0;
  23. lambda_bn[i][j] = square_front[i - 1] / volume[i][0] * 1.0 / 6.0;
  24. lambda_bf[i][j] = 0.0;
  25.  
  26. lambda_fb[i][j] = 0.0;
  27. lambda_fn[i][j] = square_front[i] / volume[i][0] * 1.0 / 6.0;
  28. lambda_ff[i][j] = 0.0;
  29.  
  30.  
  31. double prob = 1.0 - lambda_bn[i][j] - lambda_fn[i][j] - 1.0 / 6.0 * square_front[i] / volume[i][0] - 1.0 / 6.0 * square_front[i - 1] / volume[i][0];
  32. double sum = 2.0 * (square_side[i] + volume[i][0]);
  33.  
  34. lambda_nb[i][j] = prob * square_side[i] / sum;
  35. lambda_nn[i][j] = prob * 2.0 * volume[i][0] / sum;
  36. lambda_nf[i][j] = prob * square_side[i] / sum;
  37.  
  38. }
  39. }
  40.  
  41. }
  42.  
  43. void polar(vector < double > &square_front, vector < double > &square_side,
  44. vector < double > &square_up, vector <vector < double >> &volume, int Mx, double alpha, int z) { //r-layers, alpha-polar angle, h- cylindrical's height(radius)
  45. double delta_h = 1;
  46. double sum = 0, sum1 = 0;
  47.  
  48.  
  49. for (int i = 1; i < Mx + 1; ++i)
  50. square_front[i] = 0.5*alpha*(i *i - (i - 1) * (i - 1));
  51.  
  52. for (int i = 1; i < Mx + 1; ++i)
  53. square_side[i] = delta_h * (i *i - (i - 1) * (i - 1));
  54.  
  55. for (int i = 1; i < Mx + 1; ++i)
  56. square_up[i] = delta_h * i*alpha;
  57.  
  58.  
  59. for (int i = 1; i < Mx + 1; ++i) {
  60. for (int j = 1; j < alpha + 1; ++j) {
  61. volume[i][j] = 0.5*alpha*(i *i - (i - 1) * (i - 1))*delta_h;
  62. }
  63. }
  64.  
  65.  
  66. }
  67.  
  68. void transposition_torus(vector < vector < double > > &lambda_bb, vector < vector < double > > &lambda_bn,
  69. vector < vector < double > > &lambda_bf, vector < vector < double > >&lambda_fb, vector < vector < double > > &lambda_fn,
  70. vector < vector < double > > &lambda_ff, vector < vector < double > > &lambda_nb, vector < vector < double > > &lambda_nn,
  71. vector < vector < double > > &lambda_nf, vector<vector < double >> square_front, vector<vector < double >> square_left,
  72. vector<vector < double >> square_right, vector<vector < double >> square_up, vector<vector < double >> volume, double r, double sectors) {
  73.  
  74.  
  75.  
  76. for (int i = 1; i < r + 1; ++i) {
  77. for (int j = 1; j < sectors + 1; ++j) {
  78. lambda_bb[i][j] = 0;
  79. lambda_bn[i][j] = square_front[i - 1][j] / volume[i][j] * 1.0 / 6.0;
  80. lambda_bf[i][j] = 0.0;
  81.  
  82. lambda_fb[i][j] = 0.0;
  83. lambda_fn[i][j] = square_front[i][j] / volume[i][j] * 1.0 / 6.0;
  84. lambda_ff[i][j] = 0.0;
  85.  
  86.  
  87. double prob = 1.0 - 1.0 / 6.0 * square_front[i-1][j] / volume[i][j] - 1.0 / 6.0 * square_front[i - 1][j] / volume[i][j];
  88. double sum = square_left[i][j]+ square_right[i][j] + 2*volume[i][j];
  89.  
  90. lambda_nb[i][j] = prob * square_left[i][j] / sum;
  91. lambda_nn[i][j] = prob * 2.0 * volume[i][j] / sum;
  92. lambda_nf[i][j] = prob * square_right[i][j] / sum;
  93.  
  94. }
  95. }
  96. }
  97.  
  98. void torus(vector<vector < double >> &square_front, vector<vector < double >> &square_left,
  99. vector<vector < double >> &square_right, vector<vector < double >> &square_up,
  100. vector<vector < double >> &volume, double r, double sectors, double R) {
  101. double alpha = 360 / sectors;
  102. double sum1 = 0, sum2 = 0;
  103. double delta_alpha = pi * alpha / 180;
  104.  
  105.  
  106. for (int i = 1; i < r + 1; i++)
  107. for (int j = 1; j < sectors + 1; j++) {
  108. square_front[i][j] = i * delta_alpha * 2.0 * pi*(R + i * cos((2 * j - 1)*delta_alpha) / 2.0);
  109. }
  110.  
  111. for (int i = 1; i < r + 1; i++)
  112. for (int j = 1; j < sectors + 1; j++) {
  113. square_left[i][j] = (i - (i - 1)) * 2.0 * pi*(R + ((2 * i - 1) / 2.0*cos((j - 1)*delta_alpha)));
  114. }
  115.  
  116. for (int i = 1; i < r + 1; i++)
  117. for (int j = 1; j < sectors + 1; j++) {
  118. square_right[i][j] = (i - (i - 1)) * 2.0 * pi*(R + ((2 * i - 1) / 2.0 * cos(j*delta_alpha)));
  119. }
  120.  
  121. for (int i = 1; i < r + 1; i++)
  122. for (int j = 1; j < sectors + 1; j++) {
  123. square_up[i][j] = delta_alpha / 2.0*(i *i - (i - 1) * (i - 1)); //delta_alpha*(i-(i-1))*(R+(2*i-1)/2.0*cos(((2*j-1)/2.0)*delta_alpha));
  124. }
  125.  
  126. vector <double> rho(r + 1, 0);
  127.  
  128. for (int i = 1; i < r + 1; ++i) {
  129. rho[i] = 0.67*(2 * i - 1);
  130. }
  131.  
  132. for (int i = 1; i < r + 1; i++) {
  133. for (int j = 1; j < sectors + 1; j++) {
  134. volume[i][j] = square_up[i][j] * pi * 2.0*(R + rho[i] * cos(((2 * j - 1) / 2.0)*delta_alpha));
  135. }
  136. }
  137.  
  138.  
  139.  
  140. }
  141.  
  142. void spherical(double R, double r, double sectors) {
  143.  
  144. double alpha = 360 / sectors;
  145. double sum1 = 0, sum2 = 0;
  146. double delta_alpha = pi * alpha / 180;
  147.  
  148. vector < vector < double > > square_up(r + 1, vector<double>(sectors + 1, 0));
  149. vector < vector < double > > square_side(r + 1, vector<double>(sectors + 1, 0));
  150. vector < vector < double > > square_front(r + 1, vector<double>(sectors + 1, 0));
  151.  
  152. for (int i = 1; i < r + 1; i++)
  153. for (int j = 1; j < sectors + 1; j++) {
  154. square_up[i][j] = delta_alpha / 2.0*(i *i - (i - 1) * (i - 1));
  155. }
  156.  
  157. for (int i = 1; i < r + 1; i++)
  158. for (int j = 1; j < sectors + 1; j++) {
  159. square_side[i][j] = (i * i - (i - 1) * (i - 1))/2.0 * sin((j)*delta_alpha);
  160. }
  161.  
  162.  
  163. for (int i = 1; i < r + 1; i++)
  164. for (int j = 1; j < sectors + 1; j++) {
  165. square_front[i][j] = i*i * -1*cos((j-j-1)*delta_alpha);
  166. }
  167.  
  168. cout << endl << " up "<<endl;
  169.  
  170. for (int i = 1; i < r + 1; i++) {
  171. for (int j = 1; j < sectors + 1; j++)
  172. cout << square_up[i][j] << " ";
  173. cout << "\n";
  174. }
  175. cout << endl << " side " << endl;
  176.  
  177. for (int i = 1; i < r + 1; i++) {
  178. for (int j = 1; j < sectors + 1; j++)
  179. cout << square_side[i][j] << " ";
  180. cout << "\n";
  181. }
  182.  
  183. cout << endl << " front " << endl;
  184.  
  185. for (int i = 1; i < r + 1; i++) {
  186. for (int j = 1; j < sectors + 1; j++)
  187. cout << square_front[i][j] << " ";
  188. cout << "\n";
  189. }
  190.  
  191. for (int i = 1; i < sectors + 1; ++i) {
  192. sum1 += square_front[r][i];
  193. }
  194.  
  195. cout << "all square sphere " << 4.0 * pi*R*R<<" "<<sum1<<endl;
  196. cout << "all volume sphere " << 4.0 / 3.0*pi*R*R*R<<endl;
  197.  
  198.  
  199.  
  200. }
  201.  
  202. using vector3d = vector<vector<vector<double>>>;
  203.  
  204. vector < vector <double>> find_G(vector <double> u, int Mx, int My) {
  205. vector < vector < double > > G(Mx+2, vector<double>(My+2, 0));
  206.  
  207.  
  208. for (int i = 0; i < Mx+1; ++i) {
  209. for (int j = 0; j < My + 1; ++j) {
  210. G[i][j] = exp(-u[i*(My + 2) + j]);
  211. }
  212. }
  213.  
  214.  
  215. return G;
  216. }
  217.  
  218. vector3d find_Gforw(vector<vector <double>> G, vector3d &Gforw, int N,
  219. int Mx, int My, vector < vector < double > > lambda_bb, vector < vector < double > > lambda_bn,
  220. vector < vector < double > > lambda_bf, vector < vector < double > > lambda_nb,
  221. vector < vector < double > > lambda_nn, vector < vector < double > > lambda_nf,
  222. vector < vector < double > > lambda_ff, vector < vector < double > > lambda_fn,
  223. vector < vector < double > > lambda_fb) {
  224.  
  225.  
  226. for (int j = 1; j < N; ++j) {
  227. for (int k = 1; k < My + 1; ++k) {
  228. for (int i = 1; i < Mx + 1; ++i) {
  229. Gforw[i][k][j] = G[i][k] * (lambda_bb[i][k] * Gforw[i - 1][k - 1][j - 1] + lambda_bn[i][k] * Gforw[i - 1][k][j - 1] + lambda_bf[i][k] * Gforw[i - 1][k + 1][j - 1] +
  230. lambda_nb[i][k] * Gforw[i][k - 1][j - 1] + lambda_nn[i][k] * Gforw[i][k][j - 1] + lambda_nf[i][k] * Gforw[i][k + 1][j - 1] +
  231. lambda_fb[i][k] * Gforw[i + 1][k - 1][j - 1] + lambda_fn[i][k] * Gforw[i + 1][k + 1][j - 1] + lambda_ff[i][k] * Gforw[i + 1][k + 1][j - 1]);
  232. }
  233.  
  234. Gforw[Mx + 1][k][j] = Gforw[Mx][k][j]; //граничное условие
  235. Gforw[0][k][j] = 0;
  236. }
  237. }
  238.  
  239.  
  240. return Gforw;
  241. }
  242.  
  243. vector3d find_Gback(vector<vector <double>> G, vector3d &Gback, int N,
  244. int Mx, int My, vector < vector < double > > lambda_bb, vector < vector < double > > lambda_bn,
  245. vector < vector < double > > lambda_bf, vector < vector < double > > lambda_nb,
  246. vector < vector < double > > lambda_nn, vector < vector < double > > lambda_nf,
  247. vector < vector < double > > lambda_ff, vector < vector < double > > lambda_fn,
  248. vector < vector < double > > lambda_fb) {
  249.  
  250. for (int j = N - 1; j > 0; --j) {
  251. for (int k = 1; k < My + 1; ++k) {
  252. for (int i = 1; i < Mx + 1; ++i) {
  253. Gback[i][k][j - 1] = G[i][k] * (lambda_bb[i][k] * Gback[i - 1][k - 1][j] + lambda_bn[i][k] * Gback[i - 1][k][j] + lambda_bf[i][k] * Gback[i - 1][k + 1][j] +
  254. lambda_nb[i][k] * Gback[i][k - 1][j] + lambda_nn[i][k] * Gback[i][k][j] + lambda_nf[i][k] * Gback[i][k + 1][j] +
  255. lambda_fb[i][k] * Gback[i + 1][k - 1][j] + lambda_fn[i][k] * Gback[i + 1][k][j] + lambda_ff[i][k] * Gback[i + 1][k + 1][j]);
  256. }
  257.  
  258. Gback[Mx + 1][k][j - 1] = Gback[Mx + 1][N][j - 1]; //граничное условие
  259. Gback[0][k][j - 1] = 0;
  260. }
  261. }
  262.  
  263. return Gback;
  264. }
  265.  
  266. double find_q(vector3d Gforw,vector3d Gback,vector <vector <double>> volume, vector<vector < double>>G,int Mx, int My, int N) {
  267. double sum = 0, q = 0;
  268. for (int i = 1; i < Mx + 1; ++i) {
  269. for (int k = 1; k < My + 1; ++k) {
  270. sum = 0;
  271. for (int j = 0; j < N; ++j)
  272. sum += (Gforw[i][k][j] * Gback[i][k][j]) / G[i][k];
  273. q += sum*volume[i][k];
  274. }
  275. }
  276.  
  277. return q;
  278. }
  279.  
  280. vector <vector <double>> find_fi_p(vector3d Gforw, vector3d Gback, vector<vector<double>> G, double q, int Mx, int My, int N, double theta) {
  281. double sum, teta = theta * N;
  282. vector < vector < double > > fi_p(Mx + 2, vector<double>(My + 2, 0));
  283. for (int i = 1; i < Mx + 1; ++i) {
  284. for (int k = 1; k < My+1; ++k) {
  285. sum = 0;
  286. for (int j = 0; j < N; ++j)
  287. sum += (Gforw[i][k][j] * Gback[i][k][j]) / G[i][k];
  288.  
  289. fi_p[i][k] = (teta / q)*sum;
  290. }
  291. }
  292. return fi_p;
  293. }
  294. vector<vector<double>> find_fi_w(vector<vector<double>> G, int Mx, int My) {
  295. vector<vector<double>> fi_w(Mx + 1, vector<double>(My+1,0));
  296. for (int i = 1; i < Mx + 1; ++i) {
  297. for (int j = 1; j < My+1; ++j)
  298. fi_w[i][j] = G[i][j];
  299. }
  300.  
  301. return fi_w;
  302. }
  303.  
  304. vector <double> find_grad(vector<vector<double>> fi_p, vector<vector<double>> fi_w, int Mx, int My) {
  305. int M = (Mx + 2)*(My + 2);
  306. vector <double> grad(M + 1, 0);
  307. for (int i = 1; i < Mx + 1; ++i) {
  308. for (int k = 1; k < My + 1; ++k)
  309. grad[i * (My + 2) + k] = -1 + 1.0 / (fi_p[i][k] + fi_w[i][k]);
  310. }
  311.  
  312. return grad;
  313. }
  314.  
  315.  
  316. vector<double> function(vector <double> u, int My, int Mx, int N, double theta,
  317. vector<vector<double>> &fi_p, vector<vector<double>> &fi_w, double &q, double xmin,
  318. double xmax, double ymin,double ymax, vector<vector<double>> volume,
  319. vector < vector < double > > lambda_bb, vector < vector < double > > lambda_bn,
  320. vector < vector < double > > lambda_bf, vector < vector < double > > lambda_nb,
  321. vector < vector < double > > lambda_nn, vector < vector < double > > lambda_nf,
  322. vector < vector < double > > lambda_ff, vector < vector < double > > lambda_fn,
  323. vector < vector < double > > lambda_fb) {
  324.  
  325. vector3d Gback (Mx + 2, vector <vector<double>>(My + 2, vector<double>(N,0)));
  326. vector3d Gforw(Mx + 2, vector <vector<double>>(My + 2, vector<double>(N, 0)));
  327.  
  328. vector <vector <double>> G(Mx+2 , vector <double> (My+2,0));
  329. int M = (Mx + 2)*(My + 2);
  330. vector <double> grad(M, 0);
  331.  
  332. G = find_G(u, Mx, My);
  333.  
  334.  
  335. for (int i = xmin; i < xmax + 1; ++i) {
  336. for (int k = ymin; k < ymax + 1; ++k) {
  337. Gforw[i][k][0] = G[i][k];
  338. }
  339. }
  340.  
  341.  
  342. Gforw = find_Gforw(G, Gforw, N, Mx, My, lambda_bb, lambda_bn, lambda_bf, lambda_nb, lambda_nn, lambda_nf, lambda_ff, lambda_fn, lambda_fb);
  343.  
  344.  
  345.  
  346. for (int i = 1; i < Mx + 1; ++i) {
  347. for (int k = 1; k < My + 1; ++k) {
  348. Gback[i][k][N - 1] = G[i][k];
  349. }
  350. }
  351.  
  352. Gback = find_Gback(G, Gback, N, Mx, My, lambda_bb, lambda_bn, lambda_bf, lambda_nb, lambda_nn, lambda_nf, lambda_ff, lambda_fn, lambda_fb);
  353. q = find_q(Gforw, Gback, volume, G, Mx, My, N);
  354. fi_p = find_fi_p(Gforw, Gback, G, q, Mx, My, N, theta);
  355. fi_w = find_fi_w(G, Mx, My);
  356. grad = find_grad(fi_p, fi_w, Mx, My);
  357.  
  358. return grad;
  359. }
  360.  
  361. vector < double > product_matrix_on_vector(vector < vector < double > >A, vector <double> grad, int M) {
  362. vector <double> a(M + 1, 0);
  363.  
  364. for (int i = 1; i < M + 1; i++) {
  365. a[i] = 0;
  366. for (int j = 1; j < M + 1; j++) {
  367. a[i] += A[i - 1][j - 1] * grad[j];
  368. }
  369. }
  370.  
  371. return a;
  372. }
  373.  
  374.  
  375. vector < double > func_direction(vector < vector < double > > A, vector <double> grad, int M) {
  376. vector <double> direction(M + 1, 0);
  377.  
  378. direction = product_matrix_on_vector(A, grad, M);
  379.  
  380. for (int i = 0; i < M + 1; i++)
  381. direction[i] = -1 * direction[i];
  382. return direction;
  383. }
  384.  
  385. double length_of_grad(vector <double> grad, int M) {
  386. double sum = 0;
  387. for (int i = 1; i < M + 1; ++i) {
  388. sum += grad[i] * grad[i];
  389. }
  390. return sqrt(sum);
  391. }
  392.  
  393.  
  394. vector < vector < double > > product(vector < vector < double > > A, vector < vector < double > > U, int N) {
  395. vector < vector < double > > c(N, vector<double>(N, 0));
  396. for (int i = 0; i < N; i++) {
  397. for (int j = 0; j < N; j++) {
  398. c[i][j] = 0;
  399. for (int t = 0; t < N; t++)
  400. c[i][j] += A[i][t] * U[t][j];
  401. }
  402. }
  403. return c;
  404. }
  405.  
  406. vector < vector < double > > division_matrix_by_number(vector < vector < double > > matrix, double number, int N) {
  407. vector < vector < double > > matrix1(N, vector<double>(N, 0));
  408. for (int i = 0; i < N; i++) {
  409. for (int j = 0; j < N; j++)
  410. matrix1[i][j] = matrix[i][j] / number;
  411. }
  412. return matrix1;
  413. }
  414.  
  415. double make_number(vector <double> a, vector <double> b, int N) {
  416. double number = 0;
  417. for (int i = 0; i < N; i++)
  418. number += a[i] * b[i];
  419. return number;
  420. }
  421.  
  422. vector < vector < double > > make_matrix(vector <double> a, vector <double> b, int M) {
  423. vector < vector < double > > matrix(M, vector<double>(M, 0));
  424. for (int i = 0; i < M; i++) {
  425. for (int j = 0; j < M; j++)
  426. matrix[i][j] = a[i] * b[j];
  427. }
  428. return matrix;
  429. }
  430.  
  431. vector < vector < double > > formula(vector < vector < double > > matrix, vector < double > alpha, vector < double > beta, int M) {
  432. double number1 = 0, number2 = 0;
  433. vector < vector < double > > matrix1(M, vector<double>(M + 1, 0));
  434. vector < vector < double > > matrix2(M + 1, vector<double>(M + 1, 0));
  435. vector < vector < double > > B(M + 1, vector<double>(M + 1, 0));
  436. vector < vector < double > > C(M + 1, vector<double>(M + 1, 0));
  437. vector <double> a(M + 1, 0);
  438. vector <double> b(M + 1, 0);
  439. matrix1 = make_matrix(alpha, alpha, M);
  440. number1 = make_number(alpha, beta, M);
  441. B = division_matrix_by_number(matrix1, number1, M);
  442. a = product_matrix_on_vector(matrix, beta, M);
  443. matrix2 = make_matrix(a, beta, M);
  444. matrix2 = product(matrix2, matrix, M);
  445. b = product_matrix_on_vector(matrix, beta, M);
  446. number2 = make_number(b, beta, M);
  447. C = division_matrix_by_number(matrix2, number2, M);
  448.  
  449. for (int i = 0; i < M; i++) {
  450. for (int j = 0; j < M; j++)
  451. matrix[i][j] = matrix[i][j] + B[i][j] - C[i][j];
  452. }
  453.  
  454. return matrix;
  455. }
  456.  
  457.  
  458. int main() {
  459.  
  460. double theta, del, nu, max_step, xmin, xmax, ymin, ymax, q = 0, R_g = 0, k = 0;
  461. int N, Mx, My;
  462. string key;
  463.  
  464. freopen("Text.txt", "r", stdin);
  465.  
  466. cin >> N >> Mx >> My >> theta >> nu >> del >> max_step >> key >> xmin >> xmax >> ymin >> ymax;
  467.  
  468. int M = (Mx + 2)*(My + 2);
  469. double teta = N * theta;
  470.  
  471.  
  472.  
  473. vector3d Gforw(Mx + 1, vector <vector<double>>(My + 1, vector<double>(N+1, 0)));
  474. vector3d Gback(Mx + 1, vector <vector<double>>(My + 1, vector<double>(N+1, 0)));
  475. vector < vector < double > > fi_p(Mx + 1, vector<double>(My + 1,0));
  476. vector < vector < double > > fi_w(Mx + 1, vector<double>(My + 1, 0));
  477. vector <double> grad(M+2, 0);
  478. vector <double> direction(M+2, 0);
  479. vector <double> alpha(M+2, 0);
  480. vector <double> beta(M+2, 0);
  481. vector <double> u(M+2, 0);
  482. vector < vector < double > > A(M, vector<double>(M));
  483. vector < vector < double > > square_up(Mx + 2, vector<double>(My + 2, 0));
  484. vector < vector < double > > square_right(Mx + 2, vector<double>(My + 2, 0));
  485. vector < vector < double > > square_left(Mx + 2, vector<double>(My + 2, 0));
  486. vector < vector < double > > square_front(Mx + 2, vector<double>(My + 2, 0));
  487. vector < vector < double > > volume(Mx + 2, vector<double>(My + 2, 0));
  488. vector < vector < double > > lambda_bb(Mx + 2, vector<double>(My + 2, 0)); // -1, -1
  489. vector < vector < double > > lambda_bn(Mx + 2, vector<double>(My+ 2, 0)); // -1, 0
  490. vector < vector < double > > lambda_bf(Mx + 2, vector<double>(My + 2, 0)); // -1, 1
  491. vector < vector < double > > lambda_fb(Mx + 2, vector<double>(My + 2, 0)); // 1, -1
  492. vector < vector < double > > lambda_fn(Mx + 2, vector<double>(My + 2, 0)); // 1, 0
  493. vector < vector < double > > lambda_ff(Mx + 2, vector<double>(My + 2, 0)); // 1, 1
  494. vector < vector < double > > lambda_nb(Mx + 2, vector<double>(My + 2, 0)); // 0, -1
  495. vector < vector < double > > lambda_nn(Mx + 2, vector<double>(My + 2, 0)); // 0, 0
  496. vector < vector < double > > lambda_nf(Mx + 2, vector<double>(My + 2, 0)); // 0, 1
  497. vector <double> square_up1(Mx + 2);
  498. vector <double> square_front1(Mx + 2);
  499. vector <double> square_side(Mx + 2);
  500.  
  501. for (int i = 0; i < M; ++i) {
  502. for (int j = 0; j < M; ++j) {
  503. if (i == j)
  504. A[i][j] = 1;
  505. else A[i][j] = 0;
  506. }
  507. }
  508.  
  509.  
  510. if (key == "polar") {
  511.  
  512. polar(square_up1, square_side, square_front1, volume, Mx, My, N);
  513. transposition_polar(lambda_bb, lambda_bn, lambda_bf, lambda_fb, lambda_fn, lambda_ff, lambda_nb, lambda_nn,
  514. lambda_nf, square_front1, square_side, square_up1, volume, Mx, My);
  515.  
  516. }
  517. else if (key == "torus") {
  518.  
  519. torus (square_up, square_right, square_left, square_front, volume, Mx, My, N);
  520. transposition_torus(lambda_bb, lambda_bn, lambda_bf, lambda_fb, lambda_fn, lambda_ff, lambda_nb, lambda_nn,
  521. lambda_nf, square_front, square_left, square_right, square_up, volume, Mx, My);
  522.  
  523. }
  524.  
  525. grad = function(u, My, Mx, N, theta, fi_p, fi_w, q, xmin, xmax, ymin, ymax, volume, lambda_bb, lambda_bn, lambda_bf, lambda_nb,
  526. lambda_nn, lambda_nf, lambda_ff, lambda_fn, lambda_fb);
  527.  
  528. do {
  529.  
  530. direction = func_direction(A, grad, M); //to find direction
  531. /*cout << "direction ";
  532.  
  533. for (int i = 0; i < M; i++)
  534. cout << direction[i] << " ";
  535. cout << endl;*/
  536.  
  537. for (int i = 1; i < M + 1; i++) // old points
  538. alpha[i - 1] = u[i];
  539.  
  540. for (int i = 1; i < M + 1; i++) // new points
  541. u[i] = u[i] + nu * direction[i];
  542.  
  543. for (int i = 1; i < M + 1; i++) // new points - old points
  544. alpha[i - 1] = u[i] - alpha[i - 1];
  545.  
  546. for (int i = 1; i < M + 1; i++) //old values of grad
  547. beta[i - 1] = grad[i];
  548.  
  549. grad = function(u, My, Mx, N, theta, fi_p, fi_w, q, xmin, xmax, ymin, ymax, volume, lambda_bb, lambda_bn, lambda_bf, lambda_nb,
  550. lambda_nn, lambda_nf, lambda_ff, lambda_fn, lambda_fb);
  551.  
  552. for (int i = 1; i < M + 1; i++)
  553. beta[i - 1] = grad[i] - beta[i - 1];
  554.  
  555. //A = formula(A, alpha, beta, M); //Ak=Ak-1 + B - C;
  556.  
  557. k++;
  558.  
  559. cout << k << " " << length_of_grad(grad, M) << endl;
  560. } while ((length_of_grad(grad, M) > del) && (k < max_step));
  561.  
  562.  
  563.  
  564. FILE * txt = fopen("fi_p_out.txt", "w");
  565. fprintf(txt, "%5d ", (int)My);
  566. for (int j = 1; j < My + 1; j++) {
  567. fprintf(txt, " %11d", j);
  568. }
  569. fprintf(txt, "\n");
  570. for (int i = 1; i < Mx + 1; i++) {
  571. fprintf(txt, "%5d ", i);
  572. for (int j = 1; j < My + 1; j++) {
  573. fprintf(txt, " %8.5e", fi_p[i][j]);
  574. }
  575. fprintf(txt, "\n");
  576. }
  577. fclose(txt);
  578.  
  579.  
  580. _getch();
  581. return 0;
  582. }
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement
Public Pastes
Advertisement
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy.  OK, I Understand
Not a member of Pastebin yet?
Sign Up, it unlocks many cool features!