一次元イジングモデル

前に書いた２次元のモデルと比較したかったので、
一次元イジングモデルのシュミレーションを行なうコードを書いて見ました。

温度を徐々に変化させながら全体の磁化を出力させています。
結果を見ると２次元の場合に比べて、低い温度(逆温度betaが高いところ)でスピンが揃い初めることが分かります。

理論によれば 一次元の場合有限温度でスピンが揃うことは無いはずですが、系の大きさが有限であるために、beta~4程度でスピンは完全に揃ってしまいました。

#include <iostream>
#include <vector>
#include <cmath>
#include <ctime>
 
using namespace std;
 
typedef vector<int> Sv;
const int L = 10000;
const int M = 200;
const double J = 1.0;
const double B = 0.0;
double beta = 0;
Sv spin(L, 1);

double ran() {
  return (double)rand()/((double)RAND_MAX+1);
}
 
double dE(int i,int Pre) {
  return -J *( Pre - spin[i])* ( (i == 0 ? 0 : spin[i-1])
				 + (i == L-1 ? 0 : spin[i+1]))
    -B*(Pre - spin[i]);
}    
void sweep_1(int i) {
  double r = ran();
  int Pre = -spin[i];
  double difE = dE(i, Pre);
  if(difE < 0) {
    spin[i] = Pre;
  } else if(r < exp(-beta*difE) ) {
    spin[i] = Pre;
  }
}
void sweep_all() {
  for(int i = 0; i < L ; i++) {
    sweep_1(i);
  }
}
double cl(int n) {
  double sum = 0;
  for(int i = 0; i < L; i++ ) {
    sum += (double) spin[i] * spin[i];
  }
  return sum;
}
double magnet() {
  double sum = 0;
  for(int i = 0 ; i < M; i++) {
    for(int s = 0; s < L; s++) {
      sum += spin[s];
    }
    sweep_all();
  }
  return sum / (M * L);
}
void output() {
  for(beta = 0; beta < 10; beta += 0.1) {
    cout << beta << " " << magnet() << endl;
  }
}
int main(int argc, char **argv) {
  srand( (unsigned)time(NULL));
  output();
  return 0;
}