buxing.cpp

二维ixing模型的算法

#include"stdio.h"
#include"stdlib.h"
#include"math.h"
#include"time.h"

#define L 40        //矩阵尺度为L*L
#define N 1000000   //N为在温度T下的采样数    
#define N0 100000   //计算中舍去前5000个样本
#define Nt 20      //NT为温度点个数


int a[L][L]; 
long M[N];       
double E[N];
FILE *fp,*fq;

void initialize(void);
double TvsX(double B);
void randchoose(double B);
double compute(double B);

void main()
{
    double x[Nt],B;
	int i;
   //T在0.3~0.59(0.1~1.05)之间变化时，C的取值 
	for(i=0;i<Nt;i++)     
	{	
		B=0.3+i*0.01;
		printf("\nB=%f\n",B);
		
		initialize();
		x[i]=TvsX(B);
		
	}
	//将T-x数据写入文件*.txt
	fp=fopen("B.txt","w");    
	fq=fopen("x.txt","w");
   
	for(i=0;i<Nt;i++)
	{
      fprintf(fp,"%f\n",0.3+i*0.01);
	  fprintf(fq,"%f\n",x[i]);
	}

fclose(fp);
fclose(fq);
}

double TvsX(double B)
{
		
     randchoose(B);
     return compute(B);

}



void initialize(void)
{
	//矩阵初始化                         
	int i,j;
    for(i=0;i<L;i++)     
		for(j=0;j<L;j++)
           a[i][j]=1;
    //计算初始E,M  
    M[0]=0;       
	E[0]=0;    
	
    for(i=0;i<L;i++)
        for(j=0;j<L;j++)
		{
			E[0]+=-0.5*(a[i][j]*a[(L+i-1)%L][j]+a[i][j]*a[(L+i+1)%L][j]+
				      a[i][j]*a[i][(L+j-1)%L]+a[i][j]*a[i][(L+j+1)%L]);
			M[0]+=a[i][j];
		} 

}

void randchoose(double B)
{
    
	double rw;
	int rc1,rc2;
	int s;
	double E1,E2,detE;
	long i,k=1;
	double w8=exp(-8*B),w4=exp(-4*B);
	/*随机选择一个粒子改变其自旋值*/
	for(i=0;i<N-1;i++)
	 {
          rc1=rand()*(L-1)/RAND_MAX;    
          rc2=rand()*(L-1)/RAND_MAX;
		  E1=-1.0*(a[rc1][rc2]*a[(L+rc1-1)%L][rc2]+a[rc1][rc2]*a[(L+rc1+1)%L][rc2]+
			  a[rc1][rc2]*a[rc1][(L+rc2-1)%L]+a[rc1][rc2]*a[rc1][(L+rc2+1)%L]);
          s=a[rc1][rc2];  

          if(a[rc1][rc2]==1)  a[rc1][rc2]=-1;   
	       else a[rc1][rc2]=1;

		  E2=-1.0*(a[rc1][rc2]*a[(L+rc1-1)%L][rc2]+a[rc1][rc2]*a[(L+rc1+1)%L][rc2]+
			 a[rc1][rc2]*a[rc1][(L+rc2-1)%L]+a[rc1][rc2]*a[rc1][(L+rc2+1)%L]);
		  detE=E2-E1;
	     
		 bool select=false;
		 
		 if((detE==-8.0)||(detE==-4.0)) 
			 select=true;

		 if(detE==8.0)
		 {		 
			  rw=rand()*1.0/RAND_MAX;
			  if(rw<w8) 
				 select=true;
		 }
         if(detE==4.0)
		 {		 
			  rw=rand()*1.0/RAND_MAX;
			  if(rw<w4) 
				 select=true;
		 }

		 if(select==true)
		 { 
			 M[k]=M[k-1]+a[rc1][rc2]-s;
			 E[k]=E[k-1]+(E2-E1);
				 k++;
		 }
		 else 
		 {   
			 M[k]=M[k-1]; 
			 E[k]=E[k-1];
			 a[rc1][rc2]=s;
		     k++;		
		 }
	}
  

}

double compute(double B)
{
	double Ea=0,E2a=0,Ma=0,M2a=0;
	double m=0,c,x;
    long i;
	for(i=0;i<=(N-N0);i++)
	{
	
	    Ea+=E[i+N0];
	    E2a+=E[i+N0]*E[i+N0];
        Ma+=M[i+N0];
	    M2a+=M[i+N0]*M[i+N0];

	}
    
	m=Ma/(N-N0)/L/L;
    c=(E2a/(N-N0)-Ea*Ea/(N-N0)/(N-N0))*B*B;
    x=(M2a/(N-N0)-Ma*Ma/(N-N0)/(N-N0))*B;
    ////
	printf("m=%f,c=%f,x=%f\n",m,c,x);
	return x;
}