イジングモデルのシミュレーション

知り合いがやっていて楽しそうだったので僕も作ってみた。メトロポリス法というらしい。

Visual Studio 2010のWin32コンソールアプリケーション設定で作成。
1行目を消して_tmainをmainにすればg++でコンパイルできるけど。
opengl使っているのでインストールされてないとコンパイル通りません。

#include "stdafx.h"
#include <cmath>
#include <time.h>
#include <stdlib.h>
#include <GL/glut.h>

const int width=500;
const int height=500;  
const int N = 50;	// N x N boxes

const double frameratio=0.05;	// < 0.5 (枠線の幅の割合)

const double kB=1.38e-23;

const double J=1;
const double B=0;
const double Tc=J/kB*2/log(1+sqrt(2.0)); //臨界温度
const double T=3;
//const double T=Tc;


const double alpha = 0.5;
int a[N][N];
int b[N][N]; //staticに宣言しないとNが大きい時stack overflowするので

void initstates(){
	srand((unsigned int)time(NULL));
	for(int i=0;i<N;++i)for(int j=0;j<N;++j)a[i][j]=(rand()%2==0?1:-1);
}

void Initialize(){
	glutInitDisplayMode(GLUT_SINGLE | GLUT_RGB);
	glutInitWindowSize(width,height);
	glutCreateWindow("ising");
	glClearColor(1.0, 1.0, 1.0, 1.0);
	glOrtho(0, width, 0, height, -1, 1);
}

void drawfillrect_color(int x1, int y1, int x2, int y2,double r,double g,double b){
	glBegin(GL_QUADS);
		glColor3d(r,g,b);
			glVertex2d(x1,y1);
			glVertex2d(x1,y2);
			glVertex2d(x2,y2);
			glVertex2d(x2,y1);
	glEnd();
}

double E_local(int x,int y, int state){
	return -J*(a[(x+1+N)%N][y]+a[(x-1+N)%N][y]+
				a[x][(y-1+N)%N]+a[x][(y+1+N)%N])*state
			-B*state;
}

double rand1(){
	return (double)rand()/RAND_MAX;
}

void draw_boxes(){
	double xtick=(double)width/N, ytick=(double)height/N;
	double xframe=xtick*frameratio,yframe=ytick*frameratio;
	for(int i=0;i<N;++i) for (int j=0;j<N;j++){
		drawfillrect_color(
			i*xtick+xframe,j*ytick+yframe,
			(i+1)*xtick-xframe,(j+1)*ytick-yframe,
			(1+a[i][j])/2,(1-a[i][j])/2,0);
	}
}

void sweep(){
	for(int i=0;i<N;++i)for(int j=0;j<N;++j){
		if(rand1()<alpha){
			b[i][j]=a[i][j];
		}else{
			double dfE=E_local(i,j,-a[i][j])-E_local(i,j,a[i][j]);
			if(dfE<0){
				b[i][j]=-a[i][j];
			}else{
				if(rand1()<exp(-dfE/(kB*T))){
					b[i][j]=-a[i][j];
				}else{
					b[i][j]=a[i][j];
				}
			}
		}
	}
	for(int i=0;i<N;++i)for(int j=0;j<N;++j)a[i][j]=b[i][j];
}

void draw(){
	glClear(GL_COLOR_BUFFER_BIT);
	draw_boxes();
	glFlush();
}

void timer(int value){
	sweep();
	glutPostRedisplay();
	glutTimerFunc(10,timer,0);
}

int _tmain(int argc, char* argv[])
{
	glutInit(&argc,argv);
	Initialize();
	initstates();
	glutDisplayFunc(draw);
	glutTimerFunc(100,timer,0);
	glutMainLoop();
	return 0;
}

wikipediaいわくMetropolis法を発展させた最適化アルゴリズムが焼きなまし法だとかなんとか

せっかくなので大きくして枠をなくして、臨界温度で実行ししばらく放置するとこんな感じになる。特徴的なスケールのないというある意味フラクタル的な性質をもった状態。

本当ならこのあと比熱とか磁化率とかのプロットしたりして臨界指数を見積もったりするのだろうけどとりあえずここまで。