Asorander: Criou página com ' #include #include #include #include // Função para atribuir as energias do sistema double W(double w[9], double...'

2018-01-26T21:34:16Z

Criou página com '<source lang = "c"> #include <stdio.h> #include <stdlib.h> #include <math.h> #include <time.h> // Função para atribuir as energias do sistema double W(double w[9], double...'

Página nova

<source lang = "c">

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>

// Função para atribuir as energias do sistema
double W(double w[9], double);
///////////////////////////////////////////////////////////////////////////////////////////////////

// Função de atribuição de vizinhos em horizontais
int vizinhancaj(int l, int **vizinhoj)
{
int i;
for (i = 0; i < l; i++)
{
vizinhoj[i][1] = i+1;
vizinhoj[i][3] = i-1;
vizinhoj[i][2] = i;
vizinhoj[i][4] = i;
}
vizinhoj[0][3] = l-1;
vizinhoj[l-1][1] = 0;
}

// Função de atribuição de vizinhos verticais
int vizinhancai(int l, int **vizinhoi)
{
int i;
for (i = 0; i < l; i++)
{
vizinhoi[i][2] = i+1;
vizinhoi[i][4] = i-1;
vizinhoi[i][1] = i;
vizinhoi[i][3] = i;
}
vizinhoi[0][4] = l-1;
vizinhoi[l-1][2] = 0;
}

///////////////////////////////////////////////////////////////////////////////////////////////////

// Função responsável por calcular a variação de energia entre dois estados
double Eflip(int l, int** matrix, int i, int j, int **vizinhoi, int **vizinhoj)
{
int k, J = 1;
int E = 0;

for (k = 1; k < 5; k++)
{
E += matrix[vizinhoi[i][k]][vizinhoj[j][k]];
}
E = 2*J*E*matrix[i][j];
return E;
}

//////////////////////////////////////////////////////////////////////////////////////////

// Função responsável por calcular a energia total do sistema
double ETotal(int** matrix, int **vizinhoi, int **vizinhoj, int l)
{
int i, j, k, J = 1;
double E = 0;

for(i = 0; i < l; i++)
for (j = 0; j < l; j++)
{
for (k = 1; k < 5; k++)
{
E += matrix[i][j]*matrix[vizinhoi[i][k]][vizinhoj[j][k]];
}
}
E *= -J;
return E/2.0; //Divisão por 2 por somamos duas vezes cada interação entre vizinhos
}

main()
{
int i, j, k, l = 20, MCS = 0, E; // l = Lado da matrix
int t = 0;
double r, prob, M = 0, T = 2.2, ET, w[9], X, M2 = 0, M_2 = 0, C, E2, E_2; // T = Temperatura
FILE *arq, *arq2;

// Incialização das matrizes utilizadas no sistema
int ** matrix = (int **)malloc(l*sizeof(int*));
int **vizinhoi = (int **)malloc(l*sizeof(int*));
int **vizinhoj = (int **)malloc(l*sizeof(int*));

for(i=0;i<l;i++)
{
vizinhoi[i] = (int*)malloc(5*sizeof(int));
vizinhoj[i] = (int*)malloc(5*sizeof(int));
matrix[i] = (int *) malloc(l*sizeof(int*));
}

vizinhancai(l,vizinhoi);
vizinhancaj(l,vizinhoj);

srand(time(0));

///////////////////////////////////////////////////////////////////////////////////////////////////

arq = fopen("Suscetibilidade", "w");
arq2 = fopen("Calor", "w");
for(T = 2.1; T <= 2.6; T+=0.01) // Laço de variação da temperatura
{
W(w,T); // Inicializar probabilidades de mudar de estado

// Incialização da matriz de spins aleatóriamente para cada temperatura
for (i = 0; i < l; i++)
{
for (j = 0; j < l; j++)
{
r = (double)rand()/RAND_MAX;
if (r < 0.5)
//Cooperadores recebem valor 0
matrix[i][j] = (int)-1;
else
//Desertores ganham valor 1
matrix[i][j] = (int)1;
}
}

ET = ETotal(matrix, vizinhoi, vizinhoj, l); // Calculo da energia total
M = 0;

// Calculo da magnetização total
for (i = 0; i < l; i++)
for (j = 0; j < l; j++)
M += matrix[i][j];

for (MCS = 0; MCS < 1e5; MCS++) // Laço de equilíbrio do sistema
{
for (t = 0; t < l*l; t++)
{
i = rand()%l;
j = rand()%l;
E = Eflip(l, matrix, i, j, vizinhoi, vizinhoj);
if(E <= 0)
{
matrix[i][j] = -matrix[i][j];
ET += E;
M += 2*matrix[i][j];
}
else
{
prob = w[E];
r = 1.*rand()/RAND_MAX;
if(r < prob)
{
matrix[i][j] = -matrix[i][j];
ET += E;
M += 2*matrix[i][j];
}
}
}
}

// Inicialização das variáveis de médias
M2 = 0;
M_2 = 0;
E2 = 0;
E_2 = 0;

for (MCS = 0; MCS < 1e6; MCS++) // Laço de MCS
{
for (t = 0; t < l*l; t++) // Laço de tempo para 1 MCS
{

// Escolha aleatória de um spin
i = rand()%l;
j = rand()%l;

// Calculo da variação de energia
E = Eflip(l, matrix, i, j, vizinhoi, vizinhoj);

// Condicional do método de Metropolis
if(E <= 0)
{
matrix[i][j] = -matrix[i][j];
ET += E;
M += 2*matrix[i][j];
}
else
{
prob = w[E];
r = 1.*rand()/RAND_MAX;
if(r < prob)
{
matrix[i][j] = -matrix[i][j];
ET += E;
M += 2*matrix[i][j];
}
}
}

// Somatórios de média e média quadráticas por spin para magnetização e energia
M2 += (1.*M/(l*l))*(1.*M/(l*l)); // somatório da média quadrática de M
E2 += (1.*ET/(l*l))*(1.*ET/(l*l)); // somatório da média quadrática de E
E_2 += 1.*ET/(l*l); // somatório da média de E
M_2 += fabs(1.*M/(l*l)); // somatório da média de M
}

// Cálculo final das médias e médias quadráticas
M_2 = (M_2/MCS); // Média de M
E_2 = (E_2/MCS); // Média de E
M_2 *= M_2; // Quadrado da média de M
E_2 *= E_2; // Quadrado da média de E
M2 /= MCS; // Média quadrática de M
E2 /= MCS; // Média quadrática de E

X = (l*l)*(M2-M_2)/T; // Calculo da suscetibilidade magnética
C = (l*l)*(E2-E_2)/(T*T); // Calculo do calor específico

// Impressão da suscetibilidade magnética e calor específico
fprintf(arq, "%lf\t%lf\n", T, X);
fprintf(arq2, "%lf\t%lf\n", T, C);
}

fclose(arq);
fclose(arq2);
}

//Função de probabilidade de mudar de estado
double W(double w[9], double T)
{
//return exp(-1.*E/(T));
w[2] = exp(-2.0/T);
w[4] = exp(-4.0/T);
w[6] = exp(-6.0/T);
w[8] = exp(-8.0/T);
}

</source>

Suscetibilidade magnética e calor específico - Histórico de revisão

Asorander: Criou página com ' #include #include #include #include // Função para atribuir as energias do sistema double W(double w[9], double...'