SpaceQuester

Untitled

Feb 17th, 2017
343
0
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
text 3.14 KB | None | 0 0
  1. #define _USE_MATH_DEFINES
  2. #include "math.h"
  3. #include <stdlib.h>
  4. #include <stdio.h>
  5. #include <locale.h>
  6.  
  7. #define N 5;
  8. #define K 1;
  9.  
  10. const int w_rand_min = 5;
  11. const int w_rand_max = 5;
  12.  
  13. const int theta_rand_min = 0;
  14. const int theta_rand_max = 2 * M_PI;
  15.  
  16. double A[N][N];
  17. double w[N];
  18. double theta[N];
  19.  
  20. double kuramoto(int i, double theta[N])
  21. {
  22. double sum = 0;
  23.  
  24. for (int j = 0; j < N; j++)
  25. sum += A[i][j] * sin(theta[j] - theta[i]);
  26.  
  27. return w[i] + sum * (double)K / N;
  28. }
  29.  
  30. void RungeKutta(double dt, double theta[N], double theta_next[N])
  31. {
  32. double k[N][4];
  33.  
  34. // k1
  35. for (int i = 0; i < N; i++)
  36. k[i][0] = kuramoto(i, theta) * dt;
  37.  
  38. double theta_k1[N];
  39. for (int i = 0; i < N; i++)
  40. theta_k1[i] = theta[i] + k[i][0] / 2;
  41.  
  42. // k2
  43. for (int i = 0; i < N; i++)
  44. k[i][1] = kuramoto(i, theta_k1) * dt;
  45.  
  46. double theta_k2[N];
  47. for (int i = 0; i < N; i++)
  48. theta_k2[i] = theta[i] + k[i][1] / 2;
  49.  
  50. // k3
  51. for (int i = 0; i < N; i++)
  52. k[i][2] = kuramoto(i, theta_k2) * dt;
  53.  
  54. double theta_k3[N];
  55. for (int i = 0; i < N; i++)
  56. theta_k3[i] = theta[i] + k[i][2] / 2;
  57.  
  58. // k4
  59. for (int i = 0; i < N; i++)
  60. k[i][3] = kuramoto(i, theta_k3) * dt;
  61.  
  62. for (int i = 0; i < N; i++)
  63. theta_next[i] = theta[i] + (k[i][0] + 2 * k[i][1] + 2 * k[i][2] + k[i][3]) / 6;
  64. }
  65.  
  66. void CopyArray(double source[N], double target[N])
  67. {
  68. for (int i = 0; i < N; i++)
  69. target[i] = source[i];
  70. }
  71.  
  72. void main()
  73. {
  74. /*for (int i = 0; i < N; i++)
  75. {
  76. w[i] = ((double)rand() / RAND_MAX) * (w_rand_max - w_rand_min) + w_rand_min;
  77. }
  78.  
  79. for (int i = 0; i < N; i++)
  80. {
  81. theta[i] = ((double)rand() / RAND_MAX) * (theta_rand_max - theta_rand_min) + theta_rand_min;
  82. }*/
  83. w[0] = 5;
  84. w[1] = 5;
  85. w[2] = 5;
  86. w[3] = 5;
  87. w[4] = 5;
  88. theta[0] = 0;
  89. theta[1] = M_PI / 2;
  90. theta[2] = M_PI;
  91. theta[3] = (3 / 2) * M_PI;
  92. theta[4] = 2 * M_PI;
  93.  
  94.  
  95. /*{
  96. A[0][0] = 0;
  97. A[0][1] = 1;
  98. A[0][2] = 0;
  99. A[0][3] = 0;
  100. A[0][4] = 1;
  101.  
  102. A[1][0] = 1;
  103. A[1][1] = 0;
  104. A[1][2] = 1;
  105. A[1][3] = 0;
  106. A[1][4] = 0;
  107.  
  108. A[2][0] = 0;
  109. A[2][1] = 1;
  110. A[2][2] = 0;
  111. A[2][3] = 1;
  112. A[2][4] = 0;
  113.  
  114. A[3][0] = 0;
  115. A[3][1] = 0;
  116. A[3][2] = 1;
  117. A[3][3] = 0;
  118. A[3][4] = 1;
  119.  
  120. A[4][0] = 1;
  121. A[4][1] = 0;
  122. A[4][2] = 0;
  123. A[4][3] = 1;
  124. A[4][4] = 0;
  125. }*/
  126.  
  127. {
  128. A[0][0] = 1;
  129. A[0][1] = 1;
  130. A[0][2] = 1;
  131. A[0][3] = 1;
  132. A[0][4] = 1;
  133.  
  134. A[1][0] = 1;
  135. A[1][1] = 1;
  136. A[1][2] = 1;
  137. A[1][3] = 1;
  138. A[1][4] = 1;
  139.  
  140. A[2][0] = 1;
  141. A[2][1] = 1;
  142. A[2][2] = 1;
  143. A[2][3] = 1;
  144. A[2][4] = 1;
  145.  
  146. A[3][0] = 1;
  147. A[3][1] = 1;
  148. A[3][2] = 1;
  149. A[3][3] = 1;
  150. A[3][4] = 1;
  151.  
  152. A[4][0] = 1;
  153. A[4][1] = 1;
  154. A[4][2] = 1;
  155. A[4][3] = 1;
  156. A[4][4] = 1;
  157. }
  158.  
  159. const double t_max = 20;
  160. const double t_start = 0;
  161. const double dt = 0.01;
  162.  
  163. double t = t_start;
  164.  
  165. FILE *fp;
  166. fp = fopen("results.txt", "w+");
  167. //setlocale(LC_NUMERIC, "French_Canada.1252");
  168.  
  169. while (t <= t_max)
  170. {
  171. fprintf(fp, "%f\t", t);
  172. for (int i = 0; i < N; i++)
  173. {
  174. fprintf(fp, "%f\t", theta[i]);
  175. }
  176. fprintf(fp, "\n");
  177.  
  178. double theta_next[N];
  179.  
  180. RungeKutta(dt, theta, theta_next);
  181. CopyArray(theta_next, theta);
  182.  
  183. t += dt;
  184. }
  185.  
  186. fclose(fp);
  187. }
Advertisement
Add Comment
Please, Sign In to add comment