scilabバージョン - apmay @ ウィキ

//test2.sce　exec 'test2.sce'で動く（むしろ発散してすぐ黒くなる）青が正、黒が負

//パラメータ設定
N=8;//格子の大きさ
Alpha=3.0;
Beta=2.0;
%pi;
Gamma=2.0/%pi;
h_ini=1.5;
dt=0.01;//時間間隔
v=0.1;//速度
h_minus=v*dt//1ステップで減る磁場
t_0=h_ini/h_minus;//磁場が消えるステップ数


//領域確保
S=ones(N,N);//Phiを表す行列
Temp_Sa=zeros(N,N);//修正オイラー法で使用 k空間

Temp_S=zeros(N,N);//修正オイラー法で使用 k空間(追加)

Temp_Sb=zeros(N,N);//修正オイラー法で使用　k空間
Temp_SSS=zeros(N,N);//修正オイラー法で使用　k空間
T=zeros(N,N);//Phiの表示に使用する行列
SSS=zeros(N,N);//Phiの3乗にあたる行列
S_f=zeros(N,N);//Phiのフーリエ変換にあたる行列　k空間
SSS_f=zeros(N,N);//Phiの3乗のフーリエ変換にあたる行列　k空間
K=zeros(N,N);//k^2*Phi_kの計算で使う　係数
G=zeros(N,N);//G_k*Phi_kの計算で使う　係数
D=zeros(N,N);//delta_k,0にあたる行列　係数

//初期設定
t=0;//特に使いませんが時刻ということで
Step=0;//こちらも特に使いませんがStep数です
h_t=h_ini;//h(t)に相当するものです。tまたはstepごとに変化します。

  //Sの初期設定
    S=S+rand(N,N)/10;//Sの各成分が1から1.1のランダムな値をとるように設定しています。
    S_f=fft(S);//Sの高速フーリエ変換をなさってくれるありがたい組み込み関数
  //Phiの表示,Tの更新
    T=sign(real(S));//signは各成分を正ならば1、負ならば―1に変換する組み込み関数
                    //realは各成分を実部だけ取り出したものにする組み込み関数
    Matplot(10*T);//イジング模型でおなじみの行列の成分を可視化してくれるありがたい組み込み関数

//Kの定義　これはk^2にあたるものです。これを定義すると(7)式を行列計算の形で使うことができます
        //Scilabは繰り返し計算が遅いかわりに、行列計算がはやいので
for i=1:N;
  for j=1:N;
    K(i,j)=2*%pi*((i-1)^2+(j-1)^2)/N^2;
  end
end

//Gの定義　これはG_kにあたるものです
G=-2*%pi*sqrt(K)+4*ones(N,N);//G_k=4-2*pi*k

//Dの定義　δ_k,0にあたるものです。
D(1,1)=1;//他の成分はゼロになっています。


// 準備終了

//磁場がなくなるまで
for i=1:t_0;
  SSS=S.^3;
  SSS_f=fft(SSS);

  //S_fの更新　二次のルンゲクッタ法を使ったつもりです。(7)式を行列計算で更新しています。
  Temp_Sa = Alpha*S_f - Alpha*SSS_f - Beta*K.*S_f - Gamma*G.*S_f + h_t*D;
  
  Temp_S = S_f+dt*Temp_Sa;//追加
  
  Temp_SSS = ifft(Temp_S);//変更
  
  Temp_SSS = Temp_SSS.^3;
  
  Temp_SSS = fft(Temp_SSS);//この3行は、Temp_Sの逆フーリエ変換して、成分3乗して、再びフーリエ変換しました。
  
  Temp_Sb = Alpha*Temp_S - Alpha*Temp_SSS - Beta*K.*Temp_S - Gamma*G.*Temp_S + (h_t-0.5*h_minus)*D;//(修正)
  
  S_f = S_f + (Temp_Sa+Temp_Sb)*dt/2;
  
  //h(t),tの更新
  h_t=h_t-h_minus;
  t=t+dt;
  Step=Step+1;

  //Sの更新
  S=ifft(S_f);//逆フーリエ変換の組み込み関数

  
  //Phiの表示,Tの更新
  T=sign(real(S));
  Matplot(10*T);
end

//磁場がなくなった後
for i=1:t_0;
  SSS=S.^3;
  SSS_f=fft(SSS);

//S_fの更新
  Temp_Sa = Alpha*S_f - Alpha*SSS_f - Beta*K.*S_f - Gamma*G.*S_f + h_t*D;
  
  Temp_S = S_f+dt*Temp_Sa;//追加
  
  Temp_SSS = ifft(Temp_S);//変更
  
  Temp_SSS = Temp_SSS.^3;
  
  Temp_SSS = fft(Temp_SSS);//この3行は、Temp_Sの逆フーリエ変換して、成分3乗して、再びフーリエ変換しました。
  
  Temp_Sb = Alpha*Temp_S - Alpha*Temp_SSS - Beta*K.*Temp_S - Gamma*G.*Temp_S + (h_t-0.5*h_minus)*D;//(修正)
  
  S_f = S_f + (Temp_Sa+Temp_Sb)*dt/2;
  
  //tの更新
  t=t+dt;
  Step=Step+1;

  //Sの更新
  S=ifft(S_f);//逆フーリエ変換の組み込み関数


  //Phiの表示,Tの更新
  T=sign(real(S));
  Matplot(10*T);
end