scilabイジング


※上記の広告は60日以上更新のないWIKIに表示されています。更新することで広告が下部へ移動します。

Scilabで描いたイジングモデル(磁場なし)のプログラム

メインプログラム
//program Ising
//メトロポリス・アルゴリズムを用いた
//平衡状態生成器

getf 'DeltaE.sce'//あらかじめサブプログラムDelta.sceを同じディレクトリ内に作っておいてください。中身はあとに書いてあります。

L=100 //一辺に含まれるスピン(磁石)の数。好きな数字を入れてください
T=1    //温度。好きな数字を入れてください。キュリー点の理論値は2.269くらいです。
N=L^2 //スピンの合計数。


//まずランダムな初期配置を作ります
    for y = 1:L //ランダムな初期配置
      for x = 1:L
        if rand()<0.5
           spin(x,y)=1; //上向きスピン
        else
           spin(x,y)=-1; //下向きスピン
        end
      end
    end

Matplot(spin*10);//今の状態を図示してくれます。Matplotは行列の成分を可視化する組み込み関数です。このおかげで図示だけは楽です。

  //ボルツマン確率の比の計算。後で使います。
  w(1)=exp(8/T);
  w(2)=exp(4/T);
  w(3)=1;
  w(4)=exp(-4/T);
  w(5)=exp(-8/T);

  mcs=100 //モンテカルロステップ数。好きな数字あまり大きいと終わってくれません。終わらなくていやになったらScilabを閉じればいいです。



for pass =1:mcs
    //単位スピンあたりモンテカルロ・ステップの実行 これがメトロポリスアルゴリズムです。
    for ispin = 1:N
      x = int(L*rand() +1);//ここで反転させるスピン(磁石)のX座標を選びます
      y = int(L*rand() +1);//ここで反転させるスピン(磁石)のY座標を選びます。実際に反転させるかどうかは以下に出てくる乱数によります。
      dE = DeltaE(x,y);
      DD = (dE+12)/4;
      if rand()<=w(DD)    //w(DD) は反転前と反転後のexp(-E/kT)の比
        spin(x,y)=-spin(x,y);   //この操作は反転後にエネルギーが下がるならば、このスピンを必ず反転させ、そうでない場合は確率w(DD)でこのスピンを反転させる、ことを意味します。
      end
    end
  Matplot(spin*10);
  end
end


サブプログラムDeltaE

function D=DeltaE(x,y)
  //周期的境界条件あるスピンに対してその上隣り、下隣り、右隣り、左隣りを定義します。端にあるスピンについては特別な定義をします。これが周期的境界条件です。
  if x==1
    left = spin(L,y);//左端の座標にあるスピンについては左隣は存在しません。左隣りは同じ行の右端のスピンであることにします。これがこの場合の周期的境界条件です。
  else
    left = spin(x-1,y);//それ以外は一個左隣のスピンを左隣りと考えます。
  end
  if x== L
    right = spin(1,y);//同様に右端のスピンについては右隣りが存在しません。右隣りは同じ行の左端のスピンとします。
  else
    right = spin(x+1,y);//それ以外は一個右隣りのスピンが右隣りのスピンです。
  end
  if y==1
    down = spin(x,L);//下端のスピンについては、同じ列の上端のスピンを下隣りとします。
  else
    down = spin(x,y-1);
  end
  if y== L
    up = spin(x,1);//上端のスピンについては、同じ列の下端のスピンを上隣りとします。
  else
    up = spin(x,y+1);
  end
  D = 2*spin(x,y)*(left+right+down+up);//(x,y)にある
endfunction

読みにくかったかもしれませんが以上です。

3/17伊藤
上記の内容を添付ファイルにしました。
一番下にある
ising.sce
DeltaE.sce
をダウンロードしてください。
動かし方はConsoleでexec 'Ising.sce'と入力すればよいです。
添付ファイル