Ising Modelのプログラム
http://qube.phys.kindai.ac.jp/mc.pdfを参考にしてメトロポリス法でイジングモデルのシミュレーションを行ないました。
機会があれば解説記事を書きたいとお思います。とりあえずはコードだけ。
#include <iostream> #include <vector> #include <cmath> #include <ctime> using namespace std; typedef vector<int> Sv; typedef vector<Sv> Sm; const int L = 20; const int M = 200; const double J = 1.0; const double B = 0.0; double beta = 0.01; Sm spin(L, Sv(L, 1)); double ran() { return (double)rand()/((double)RAND_MAX+1); } double E(Sm spin) { double sum = 0.0; for(int i = 0; i < L-1 ; i++) { for( int j = 0; j < L-1 ; j++) { sum += -J*(spin[i][j]*spin[i][j+1] + spin[i][j]*spin[i+1][j]); } } for(int i = 0; i < L; i++) { for(int j = 0; j < L ; j++) { sum += -B*spin[i][j]; } } return sum; } double dE(int i, int j, int Pre) { return -J *( Pre - spin[i][j] )* ( (i == 0 ? 0 : spin[i-1][j]) + (i == L-1 ? 0 : spin[i+1][j]) + (j == L-1 ? 0 : spin[i][j+1]) + (j == 0 ? 0 : spin[i][j-1])) -B*(Pre - spin[i][j]); } void sweep_1(int i, int j) { double r = ran(); int Pre = -spin[i][j]; double difE = dE(i, j, Pre); if(difE < 0) { spin[i][j] = Pre; } else if(r < exp(-beta*difE) ) { spin[i][j] = Pre; } } void sweep_all() { for(int i = 0; i < L ; i++) { for(int j = 0; j < L ; j++) { sweep_1(i, j); } } } void output() { for(beta = 0; beta < 1; beta += 0.01) { int sum = 0; for(int i = 0; i < L; i++ ) { for(int j = 0; j < L; j++) { spin[i][j] = pow(-1.0, rand()); } } for(int i = 1; i < M ; i++) { for(int j = 0; j< L; j++) { for(int k = 0; k < L; k++) { sum += spin[j][k]; } } sweep_all(); } cout << sum / M << endl; } } int main() { srand( (unsigned)time(NULL)); output(); return 0; }