「scilabイジング」の編集履歴（バックアップ）一覧はこちら
「scilabイジング」(2009/03/17 (火) 15:27:33) の最新版変更点
追加された行は緑色になります。
削除された行は赤色になります。
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伊藤
上記の内容を添付ファイルにしました。
&ref(ising.sce)
&ref(DeltaE.sce)
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'と入力すればよいです。
