C语言二维伊辛模型的仿真

话不多说，上代码

#include<stdio.h>
#include<math.h>
#include<stdlib.h>
#include<time.h>
int main()
{
//初始相位 
    float J = 1.0; //耦合常数J（勿动！） 
    int N = 20; //伊辛模型的线度
    float T,delta,R ; //模型温度T，取值0~5 ,接受率为R
    float KB = 1.0; // 玻尔兹曼常量 
	int a[N][N];//模型 
	FILE *fp = NULL;	
	FILE *fz = NULL;
	fp=fopen("磁化强度&温度.txt","a");	
	fz=fopen("能量&时间.txt","a");
	
	srand((unsigned)time(NULL));
	int i,j,xun;
	float EX=0;
	float MX=0;
	float E=0;

	int bu=5000;
	float Mx=0;
	float MM=0;
	float EE=0;
	
	for(T=0.1;T<=4.5;T+=0.01){ 
	    for(i=0;i<N;i++){
		for(j=0;j<N;j++){
			a[i][j] = rand()%2;
			if(a[i][j]==0){
				a[i][j] = -1 ;
			}}}
	    for(xun=1;xun<=bu;xun+=1){
	      for(i=0;i<N;i++){
		    for(j=0;j<N;j++){
		    	E = (a[i][(abs(j+19))%20] + a[i][(abs(j+1))%20] + a[(abs(i+19))%20][j] + a[(abs(i+1))%20][j])*(-J)*a[i][j];	
			a[i][j]= -a[i][j];
			    EX = (a[i][(abs(j+19))%20] + a[i][(abs(j+1))%20] + a[(abs(i+19))%20][j] + a[(abs(i+1))%20][j])*(-J)*a[i][j];	
			a[i][j]= -a[i][j];//返回最初值，待判断。 
			delta = E - EX;//能量差值
			R = exp(delta/(KB*T));//偏转率 
			float suiji;
			suiji = rand()*1.0/RAND_MAX;//产生随机数
			if(R>suiji||R>=1){
				a[i][j] = -a[i][j];}}
		}}
		for(xun=1;xun<=1000;xun+=1){
	      for(i=0;i<N;i++){
		    for(j=0;j<N;j++){
		    	E = (a[i][(abs(j+19))%20] + a[i][(abs(j+1))%20] + a[(abs(i+19))%20][j] + a[(abs(i+1))%20][j])*(-J)*a[i][j];	
			a[i][j]= -a[i][j];
			    EX = (a[i][(abs(j+19))%20] + a[i][(abs(j+1))%20] + a[(abs(i+19))%20][j] + a[(abs(i+1))%20][j])*(-J)*a[i][j];	
			a[i][j]= -a[i][j];//返回最初值，待判断。 
			delta = E - EX;//能量差值
			R = exp(delta/(KB*T));//偏转率 
			float suiji;
			suiji = rand()*1.0/RAND_MAX;//产生随机数
			if(R>suiji||R>=1){
				a[i][j] = -a[i][j];}
		}}
		Mx=0;
		for(i=0;i<N;i++){
		for(j=0;j<N;j++){
		  Mx += (float)a[i][j]/(float)(N*N);	
		  }}
		  MM += fabs(Mx);} 
	printf("\n");
    printf("平衡态磁化强度平均值为：%f,温度为：%f",MM/1000,T);
    fprintf(fp,"%f",MM/1000);
    fputc('\n',fp);
    MM=0;
    E=0;
    EE=0;
}
fclose(fp);
float Ex;
//磁化强度随时间的变化 
T = 2.0;
for(i=0;i<N;i++){
		for(j=0;j<N;j++){
			a[i][j] = rand()%2;
			if(a[i][j]==0){
				a[i][j] = -1 ;
			}}}
for(xun=1;xun<=1500;xun++){
        for(i=0;i<N;i++){
		    for(j=0;j<N;j++){
		    	E = (a[i][(abs(j+19))%20] + a[i][(abs(j+1))%20] + a[(abs(i+19))%20][j] + a[(abs(i+1))%20][j])*(-J)*a[i][j];	
			a[i][j]= -a[i][j];
			    EX = (a[i][(abs(j+19))%20] + a[i][(abs(j+1))%20] + a[(abs(i+19))%20][j] + a[(abs(i+1))%20][j])*(-J)*a[i][j];	
			a[i][j]= -a[i][j];//返回最初值，待判断。 
			delta = E - EX;//能量差值
			R = exp(delta/(KB*T));//偏转率 
			float suiji;
			suiji = rand()*1.0/RAND_MAX;//产生随机数
			if(R>suiji||R>=1){
				a[i][j] = -a[i][j];}}}	
        Ex=0;
		for(i=0;i<N;i++){
		for(j=0;j<N;j++){
		  Ex += (a[i][(abs(j+19))%20] + a[i][(abs(j+1))%20] + a[(abs(i+19))%20][j] + a[(abs(i+1))%20][j])*(-J)*a[i][j];}}
		printf("循环次数:%d,能量：%f\n",xun,Ex);
		fprintf(fz,"%f",Ex);
		fputc('\n',fz);	
	}
fclose(fz);
} 

欢迎评论交流。