Advertisement
Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- #define _USE_MATH_DEFINES
- #include "math.h"
- #include <stdlib.h>
- #include <stdio.h>
- #include <locale.h>
- #include <time.h>
- #define N 10
- #define K 2.9
- const int omega_rand_glial_min = 0;
- const int omega_rand_glial_max = 1;
- const int omega_rand_neuron_min = 5;
- const int omega_rand_neuron_max = 6;
- const int theta_rand_min = 0;
- const int theta_rand_max = 2 * M_PI;
- double A[N][N];
- double omega[N];
- double theta[N];
- double kuramoto(int i, double theta[N])
- {
- double sum = 0;
- for (int j = 0; j < N; j++)
- sum += A[i][j] * sin(theta[j] - theta[i]);
- return omega[i] + sum * (double)K / N;
- }
- void RungeKutta(double dt, double theta[N], double theta_next[N])
- {
- double k[N][4];
- // k1
- for (int i = 0; i < N; i++)
- k[i][0] = kuramoto(i, theta) * dt;
- double theta_k1[N];
- for (int i = 0; i < N; i++)
- theta_k1[i] = theta[i] + k[i][0] / 2;
- // k2
- for (int i = 0; i < N; i++)
- k[i][1] = kuramoto(i, theta_k1) * dt;
- double theta_k2[N];
- for (int i = 0; i < N; i++)
- theta_k2[i] = theta[i] + k[i][1] / 2;
- // k3
- for (int i = 0; i < N; i++)
- k[i][2] = kuramoto(i, theta_k2) * dt;
- double theta_k3[N];
- for (int i = 0; i < N; i++)
- theta_k3[i] = theta[i] + k[i][2] / 2;
- // k4
- for (int i = 0; i < N; i++)
- k[i][3] = kuramoto(i, theta_k3) * dt;
- for (int i = 0; i < N; i++)
- theta_next[i] = theta[i] + (k[i][0] + 2 * k[i][1] + 2 * k[i][2] + k[i][3]) / 6;
- }
- void CopyArray(double source[N], double target[N])
- {
- for (int i = 0; i < N; i++)
- target[i] = source[i];
- }
- int main()
- {
- srand(time(NULL));
- for (int i = 0; i < N/2; i++)
- {
- omega[i] = ((double)rand() / RAND_MAX) * (omega_rand_glial_max - omega_rand_glial_min) + omega_rand_glial_min;
- }
- for (int i = N/2; i < N; i++)
- {
- omega[i] = ((double)rand() / RAND_MAX) * (omega_rand_neuron_max - omega_rand_neuron_min) + omega_rand_neuron_min;
- }
- /*w[0] = 5;
- w[1] = 5;
- w[2] = 5;
- w[3] = 5;
- w[4] = 5;*/
- for (int i = 0; i < N; i++)
- {
- theta[i] = ((double)rand() / RAND_MAX) * (theta_rand_max - theta_rand_min) + theta_rand_min;
- }
- /*theta[0] = 0;
- theta[1] = M_PI / 2;
- theta[2] = M_PI;
- theta[3] = ((double)3 / 2) * M_PI;
- theta[4] = 2 * M_PI;*/
- /*{
- A[0][0] = 0;
- A[0][1] = 1;
- A[0][2] = 0;
- A[0][3] = 0;
- A[0][4] = 1;
- A[1][0] = 1;
- A[1][1] = 0;
- A[1][2] = 1;
- A[1][3] = 0;
- A[1][4] = 0;
- A[2][0] = 0;
- A[2][1] = 1;
- A[2][2] = 0;
- A[2][3] = 1;
- A[2][4] = 0;
- A[3][0] = 0;
- A[3][1] = 0;
- A[3][2] = 1;
- A[3][3] = 0;
- A[3][4] = 1;
- A[4][0] = 1;
- A[4][1] = 0;
- A[4][2] = 0;
- A[4][3] = 1;
- A[4][4] = 0;
- }*/
- /*{
- A[0][0] = 1;
- A[0][1] = 1;
- A[0][2] = 1;
- A[0][3] = 1;
- A[0][4] = 1;
- A[1][0] = 1;
- A[1][1] = 1;
- A[1][2] = 1;
- A[1][3] = 1;
- A[1][4] = 1;
- A[2][0] = 1;
- A[2][1] = 1;
- A[2][2] = 1;
- A[2][3] = 1;
- A[2][4] = 1;
- A[3][0] = 1;
- A[3][1] = 1;
- A[3][2] = 1;
- A[3][3] = 1;
- A[3][4] = 1;
- A[4][0] = 1;
- A[4][1] = 1;
- A[4][2] = 1;
- A[4][3] = 1;
- A[4][4] = 1;
- }*/
- FILE *fp0;
- fp0 = fopen("A.txt", "w+");
- for (int i = 0; i < N; i++)
- {
- for (int j = 0; j < N; j++)
- {
- A[i][j] = 1;
- /*if (i == j)
- {
- A[i][j] = 0;
- }*/
- printf("i: %d\n", i);
- printf("j: %d\n", j);
- printf("A: %d\n", A[i][j]);
- fprintf(fp0, "%d\t", A[i][j]);
- }
- fprintf(fp0, "\n");
- }
- fclose(fp0);
- const double t_start = 0;
- const double t_max = 25;
- const double dt = 0.1;
- double t = t_start;
- FILE *fp1;
- fp1 = fopen("omega.txt", "w+");
- for (int i = 0; i < N; i++)
- {
- fprintf(fp1, "%f\n", omega[i]);
- }
- fclose(fp1);
- FILE *fp2;
- fp2 = fopen("results.txt", "w+");
- //setlocale(LC_NUMERIC, "French_Canada.1252");
- while (t < t_max + dt)
- {
- fprintf(fp2, "%f\t", t);
- for (int i = 0; i < N; i++)
- {
- fprintf(fp2, i == N - 1 ? "%f" : "%f\t", theta[i]);
- }
- fprintf(fp2, "\n");
- double theta_next[N];
- RungeKutta(dt, theta, theta_next);
- CopyArray(theta_next, theta);
- t += dt;
- printf("Progress: %d%%\n", (int)(100 * t / t_max));
- }
- fclose(fp2);
- }
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement