Scilabでイジング模型 その4

Scilabでイジング模型 その2では一次元のイジング模型の磁区が形成されていく様子をシミュレーションしました。
今回はイジング模型のモンテカルロシミュレーションから巨視的な物理量(エネルギーや磁化など)の計算を行います。
巨視的な系のエネルギー<E>は、状態|αj>におけるエネルギーE(αj)の統計平均として以下のように表されます。

\langle E \rangle = \sum_{\alpha_j}E(\alpha_j)P(\alpha_j)
P(\alpha_{j})=\frac{1}{Z(T)}\exp\left(- \frac{E(\alpha_{j})}{kT} \right)
Z(T)=\sum_{\alpha_j}\exp\left(-\frac{E(\alpha_j)}{kT}\right)

しかしながら、実際には分配関数Z(T)を計算するのが困難であるので、P(αj)に比例した確率で|αj>が出現するようなモンテカルロシミュレーションをm回行い、その平均から

\langle E \rangle = \frac{1}{m} \sum_{t=1}^{m}E_{t}(\alpha)

のように計算を行います。


一次元イジング模型における巨視的な物理量


Scilabでイジング模型 その2では一次元のイジング模型の磁区が形成されていく様子をシミュレーションしました。
今回はイジング模型のモンテカルロシミュレーションから巨視的な物理量(エネルギーや磁化など)の計算を行います。

n粒子スピン系におけるある状態ベクトル|αj>における系のエネルギーは以下のようになることをScilabでイジング模型 その2にも書きました。

E(\alpha_{j}) = - J \sum_{i=1}^{n-1}s_{i}s_{i+1}

よってある温度Tのときの状態|αj>が決まれば、系のエネルギーが決まることになりますが、実際には一意に決まるわけではなく、ボルツマン因子exp(-E(αj)/kT))に比例した確率P(αj)で色々な状態を取り得ます。

P(\alpha_{j})=\frac{1}{Z(T)}\exp\left(- \frac{E(\alpha_{j})}{kT} \right)
Z(T)=\sum_{\alpha_j}\exp\left(-\frac{E(\alpha_j)}{kT}\right)

ここでZ(T)は分配関数と呼ばれすべての状態の和です。

ひとたび確率P(αj)が求まれば、巨視的な物理量<A>は状態|αj>における物理量A(αj)を用いて以下のようにあらわすことができます。

\langle A \rangle = \sum_{\alpha_j}A(\alpha_j)P(\alpha_j)

例えばエネルギーEの場合は以下のようになります。

\langle E \rangle = \sum_{\alpha_j}E(\alpha_j)P(\alpha_j)

しかし実際にZ(T)をすべて計算するのは不可能です。代わりにP(αj)に比例した確率で|αj>を出現させるアルゴリズムがメトロポリスのアルゴリズムでした。

そこで系の巨視的なエネルギー(等の物理量)を計算する際には、メトロポリスのアルゴリズムを複数回繰り返して得られた状態に対するエネルギーの平均を代わりに用います。m回繰り返す場合は

\langle E \rangle = \frac{1}{m} \sum_{t=1}^{m}E_{t}(\alpha)

となります。

他にも微分方程式による物理現象のモデル化(PDF)には磁化M、比熱C、磁化率χの表式が以下のように与えられています。

M(\alpha_j)=\sum_{i=1}^{n}s_i
\langle M \rangle = \frac{1}{m} \sum_{t=1}^{m}M_{t}(\alpha)

C = \frac{\langle E^2 \rangle - \langle E \rangle^2}{kT}

\chi = \frac{\langle M^2 \rangle - \langle M \rangle^2}{kT}

Scilabスクリプト


これらを踏まえたScilabスクリプトを以下に示します。

clear;

// *** 定数の設定 ***
n = 100; // 粒子の数
m = 200; // 熱力学的な平均を取る回数
J = 1; // 交換エネルギー
rand("uniform"); // 乱数は一様乱数とする
tmax = 5 * n; // 時間の最大ステップ
h = 0.0; // 外部磁場
// 温度
ktmin = 0.5; // 最低温度
ktmax = 5.0; // 最高温度
nkt = 19; // 温度の分割数
//T = linspace(ktmin, ktmax, nkt); // 低温から開始
//spin = ones(1,n); // 各粒子におけるスピン(コールドスタート)
T = linspace(ktmax, ktmin, nkt); // 高温から開始
spin = 1 - 2 * round(rand(1,n)); // 各粒子におけるスピン(ランダム)

// *** エネルギーの計算関数 ***
function e = energy(spin)
e = - J * sum(spin .* [spin(2:n), spin(1)]) - h * sum(spin);
endfunction

// *** 行列の初期化 ***
E = []; // エネルギーの和
E2 = []; // エネルギーの二乗の和
M = []; // 磁化の和
M2 = []; // 磁化の二乗の和

// *** 温度のループ ***
for kt = 1:nkt do
// エネルギーの初期化
ene1 = 0; // エネルギーの和
ene2 = 0; // エネルギーの二乗和
// 磁化の初期化
mag1 = 0; // 磁化の和
mag2 = 0; // 磁化の二乗和f
// *** 熱力学平均のループ ***
for samp = 1:m do
// *** 時間発展のループ ***
for t = 1:tmax do
oldenergy = energy(spin);
element = ceil(n * rand()); // 粒子を一つ選ぶ
spin(element) = -1 * spin(element); // スピンを反転
newenergy = energy(spin);
if (newenergy > oldenergy) & (exp((- newenergy + oldenergy) / T(kt)) < rand()) then
spin(element) = -1 * spin(element); // 棄却
end
end
ene1 = ene1 + energy(spin); // エネルギーの和
ene2 = ene2 + energy(spin)^2; // エネルギーの二乗の和
mag1 = mag1 + sum(sum(spin)); // 磁化の和
mag2 = mag2 + sum(sum(spin))^2; // 磁化の二乗和
end
E = [E, ene1 / m]; // エネルギーの和
E2 = [E2, ene2 / m]; // エネルギーの二乗の和
M = [M, mag1 / m]; // 磁化の和
M2 = [M2, mag2 / m]; // 磁化の二乗和
end

// *** エネルギーと磁化の揺らぎ ***
C = (E2 - E .^ 2) ./ (n * T .^ 2); // 比熱
X = (M2 - M .^ 2) ./ (n * T); // 磁化率

// *** 厳密解の計算 ***
// 温度ベクトル
Ta = linspace(0.1,5,50);
// 粒子1個あたりの平均エネルギー
Ea = - tanh(J ./ Ta);
// 比熱
Ca = (J ./ Ta) .^ 2 ./ cosh(J ./ Ta) .^ 2;
// 磁化
Ma = sinh(h ./ Ta) ./ sqrt(sinh(h ./ Ta) .^ 2 + exp(-4 * J ./ Ta));
// 磁化率
Xa = exp(2 * J ./ Ta) ./ Ta;
RXa = 1 ./ Xa;

// *** グラフのプロット ***
// エネルギー
subplot(2,2,1);
plot(T, E ./ n, 'or');
plot(Ta, Ea, '--g');
xlabel("kT/J");
ylabel("E/NJ");
// 比熱
subplot(2,2,2);
plot(T, C, 'or');
plot(Ta,Ca,'--g');
xlabel("kT/J");
ylabel("C/Nk");
// 磁化
subplot(2,2,3);
plot(T, M ./ n, 'or');
plot(Ta,Ma ./ n,'--g');
xlabel("kT/J");
ylabel("M/N");
// 磁化率
subplot(2,2,4);
plot(T, 1 ./ X, 'or');
plot(Ta,RXa,'--g');
xlabel("kT/J");
ylabel("N/JX");


結果はScilabでイジング模型 その1と同様になります(なので下記の画像は使いまわしです)。コードもほとんど同じです。

001_20141130214554d90.png

Fig.1: 一次元のイジングモデル


関連エントリ




参考URL




付録


このエントリで使用したScilabのシミュレーション用ファイルを添付します。ファイル名末尾の".txt"を削除して、"_"を"."に変更すれば使えるはずです。(参考:ねがてぃぶろぐの付録)


参考文献/使用機器




フィードバック



にほんブログ村 その他趣味ブログ 電子工作へ

 ↑ 電子工作ブログランキング参加中です。1クリックお願いします。


コメント・トラックバックも歓迎です。 ↓      


 ↓ この記事が面白かった方は「拍手」をお願いします。

tag: Scilab 確率論 乱数 イジング模型 モンテカルロ解析 強磁性 反強磁性 

Scilabでイジング模型 その3

Scilabでイジング模型 その2では一次元のイジング模型を用いて磁区が形成されていく様子を確認しました。今回のエントリでは、二次元のイジングモデルで同様の計算をおこないます。

001_20150125120301efd.gif

Fig.1: 二次元イジング模型の時間発展。ランダムな初期状態から次第に磁区が形成されていく様子が観察できる。



二次元のイジング模型


Scilabでイジング模型 その2では計算物理学 応用編ising.cをScilabへ移植しました。今回は、このプログラムを二次元へ拡張します。

二次元の場合、スピンの向きも二次元にする(XY模型)ことも可能ですが、今回は一次元の場合と同様にスピンの向きは二つの状態しかとらないイジング模型として扱うことにします。更に相互作用する粒子は、隣接するものだけであるという仮定もそのまま使うと、主な変更点はスピンを保存する変数であるspinをベクトルから二次元の行列へと変更するところという事になります。
これに伴ってエネルギーを計算する関数も二次元へと拡張します。具体的には、列方向と行方向の両方に対してシフトと乗算を行うというだけですが。

Scilabスクリプト


二次元のイジングモデルのScilabスクリプトはising2d_sce.txtとなりました。

clear;

// *** 定数の設定 ***
n = 20; // 粒子の数
kt = 1.0; // 温度
J = 1; // 交換エネルギー (1: 強磁性, -1:反強磁性)
rand("uniform"); // 乱数は一様乱数とする
tmax = 5000; // 時間の最大ステップ

// *** 初期化 ***
// 各粒子におけるスピン
//spin = ones(n,n); // コールドスタート
spin = 1 - 2 * round(rand(n,n));

isoview(0,n,0,n)

// *** エネルギーの計算関数 ***
function e = energy(spin)
e = - J * sum(spin .* [spin(:,2:n), spin(:,1)]) - J * sum(spin .* [spin(2:n,:); spin(1,:)]);
endfunction

// *** 時間発展 ***
for t = 1:tmax do
oldenergy = energy(spin);
// 粒子を一つ選ぶ
elementx = ceil(n * rand());
elementy = ceil(n * rand());
// スピンを反転
spin(elementx,elementy) = -1 * spin(elementx,elementy);
newenergy = energy(spin);
spin(elementx, elementy) = (- 2 * ((newenergy > oldenergy) & (exp((- newenergy + oldenergy) / kt) <= rand())) + 1) * spin(elementx, elementy);
// スピン状態をプロット
Matplot((spin + 1) .* 15);
end


磁区の形成


強磁性的な相互作用を持つ系に対して、ランダムなスピンをもつ初期状態からスタートすると、終状態までに磁区が形成されていく様子を観察することができます。

002_2015012512030057b.png

003_20150125120300b77.png
Fig.2-3: 初期状態と終状態。ランダムなスピン分布をもっていた初期状態から、磁区が形成された終状態まで時間発展する。


gif動画のためのスクリプト


各時間発展サイクルで、ディスプレイにグラフを描きだす代わりに、画像出力を行うようにしておき、あとからgimpなどの画像編集ソフトを用いることでgif動画を作成することができます。ising2d-gif_sce.txt

関連エントリ




付録


このエントリで使用したScilabのシミュレーション用ファイルを添付します。ファイル名末尾の".txt"を削除して、"_"を"."に変更すれば使えるはずです。(参考:ねがてぃぶろぐの付録)


参考文献/使用機器



フィードバック



にほんブログ村 その他趣味ブログ 電子工作へ

 ↑ 電子工作ブログランキング参加中です。1クリックお願いします。


コメント・トラックバックも歓迎です。 ↓      


 ↓ この記事が面白かった方は「拍手」をお願いします。

tag: Scilab 確率論 乱数 イジング模型 モンテカルロ解析 強磁性 反強磁性 

Scilabで床で跳ねるボール

matlab 衝突のシミュレーションで質問された床ではねるボールのシミュレーションをScilabで行いました。

001_2015012510244301d.png
Fig.1: 跳ねるボールのシミュレーション


この問題は、常微分方程式ソルバodeと非線形方程式ソルバfsolveの二つを組み合わせることによってシミュレーションすることができます。


問題設定


高さx0から初速度v0でボールを投げる。ボールは重力加速度gで下方へ引かれ、床に落下する。落下したボールは、反発係数kで上方へ跳ねるとする。

この問題はmatlab 衝突のシミュレーション跳ねるボールのシミュレーションです。今回のエントリでは、この問題をScilabでシミュレーションします。

常微分方程式ソルバ


まず前半部分の放物運動は、常微分方程式ソルバodeを用いて簡単に計算できます。(参考: 常微分方程式タグ)

解くべき連立常微分方程式は以下のようになります。

\frac{\mathrm{d}x}{\mathrm{d}t} = v
\frac{\mathrm{d}v}{\mathrm{d}t} = -g

非線形方程式ソルバ


次に床で跳ねる計算を行います。
床ではねる直前の速度をv-,直後の速度をv+とすると反発係数kを用いて

v+ = - k v-

となります。

前述の常微分方程式を非線形方程式ソルバfsolveを利用して x=0 となる時刻を計算し、それぞれの跳ねるイベントの間の計算を行えば、この問題を解くことができます。

プログラム


結局Scilabスクリプトはball_sce.txtとなりました。

clear;

// *** 計算条件 ***
g = 9.81; // 重力加速度
k = 0.8; // 反発係数
// 時間
ts = 0; // 開始時刻
te = 10; // 終了時刻
// 初期条件
x0 = 10; // 初期位置
v0 = 15; // 初速度

// *** 常微分方程式の定義 ***
function dx = fall(t,x)
dx(1) = x(2); // dx/dt = v
dx(2) = -g; // dv/dt = -g
endfunction

// *** 解くべき非線形方程式 ***
function x = bound(t)
X = ode([x0; v0], ts, t, fall);
x = X(1);
endfunction

// *** メイン ***
// 初期化
tb = ts; // 衝突時刻
X = []; // プロット用の位置と速度
T = []; // プロット用の時間
// 終了時刻まで繰り返し
while tb < te
// 衝突時刻を計算
tb = fsolve(2 * te, bound);
// プロット用の時間ベクトルを作成
if tb > te then
// 終了時刻まで
time = linspace(ts, te);
else
// 衝突時刻まで
time = linspace(ts, tb);
end
T = [T, time];
// プロット用の位置と速度を計算
X = [X, ode([x0; v0], ts, time, fall)];
// 次の繰返しのための初期条件
x0 = 0; // 初期位置
v0 = -k * X(2,$); // 初速度
ts = tb; // 計算開始時刻
end

// *** グラフのプロット ***
// 位置のプロット
subplot(2,1,1);
plot(T, X(1,:),'g');
xgrid(color("gray"));
xlabel("Time");
ylabel("Position");
// 速度のプロット
subplot(2,1,2);
plot(T, X(2,:),'r');
xgrid(color("gray"));
xlabel("Time");
ylabel("Velocity");


関連エントリ




参考URL




付録


このエントリで使用したScilabのシミュレーション用ファイルを添付します。ファイル名末尾の".txt"を削除して、"_"を"."に変更すれば使えるはずです。(参考:ねがてぃぶろぐの付録)


参考文献/使用機器




フィードバック



にほんブログ村 その他趣味ブログ 電子工作へ

 ↑ 電子工作ブログランキング参加中です。1クリックお願いします。


コメント・トラックバックも歓迎です。 ↓      


 ↓ この記事が面白かった方は「拍手」をお願いします。

tag: Scilab 非線形方程式ソルバ 常微分方程式 

AkaiKKRの実行コマンド

AkaiKKR(machikaneyama)の実行ファイルspecxは通常のコマンドラインプログラムなので、リダイレクトやパイプを利用していろいろ便利なことができます。
今回のエントリでは、出力結果をファイルへリダイレクトする方法やその内容をリアルタイムに確認する方法についてまとめました。


AkaiKKRの実行


AkaiKKR(machikaneyama)の実行ファイルはspecxです。その実行は、通常のコマンドラインプログラムと同様にリダイレクトを利用します。

specx < in/infile

上記のようにすれば、結果が標準出力に表示されます。

specx < in/infile > out/outfile & 

上記のようにすれば、結果がout/outfileにリダイレクトされます。

specx < in/infile >> out/outfile & 

上記のようにリダイレクト記号を >> にすれば、結果を上書きせずにファイルの末尾に追加するようになります。

specxを実行している最中であっても下記のようにcatコマンドを使って出力ファイルの中身を確認できます。

cat out/outfile


tailコマンドで監視


さてデータが書き込まれるごとにcatコマンドで確認するのもよいですが、やはりリアルタイムに出力結果を確認したいという要求もあります。そこでtailコマンドを利用します。

specx < in/infile > out/outfile &
tail -f out/outfile

tailコマンドは、ファイルの末尾を表示するコマンドですが、-fオプションを付けることでファイルの内容を監視することができます。tailコマンドを終了するにはCtrl-Cを押します。

ただし、環境によっては出力ファイルへの書き込みがリアルタイムに行われず、ある程度バッファしてから行われることがあるようです。そのような環境では次に示すteeコマンドを利用します。

teeコマンドでファイルとプロンプトの両方に出力


teeコマンドは、標準入力から読み込んだ内容を、標準出力とファイルの両方に出力します。specxの出力はパイプを利用して渡します。

specx < in/infile | tee out/outfile

なお、出力ファイルを上書きせずに末尾に追加するには-aオプションを利用します。

specx < in/infile | tee -a out/outfile


こちらもCtrl-Cで終了しますが、そのばあいspecxまでついでに終了してしまいます。

関連エントリ




参考URL




参考文献/使用機器




フィードバック



にほんブログ村 その他趣味ブログ 電子工作へ

 ↑ 電子工作ブログランキング参加中です。1クリックお願いします。


コメント・トラックバックも歓迎です。 ↓      


 ↓ この記事が面白かった方は「拍手」をお願いします。

tag: AkaiKKR machikaneyama シェルスクリプト 

FC2カウンター
カテゴリ
ユーザータグ

LTspiceAkaiKKRmachikaneyamaScilabKKRPSoCOPアンプCPA強磁性PICモンテカルロ解析常微分方程式odeトランジスタecalj状態密度DOSインターフェース定電流スイッチング回路PDS5022半導体シェルスクリプト乱数レベルシフトHP6632A温度解析分散関係I2Cトランジスタ技術R6452A可変抵抗ブレッドボードセミナーバンドギャップ数値積分確率論反強磁性偏微分方程式バンド構造絶縁熱設計非線形方程式ソルバフォトカプラシュミットトリガLEDLM358カオスISO-I2C三端子レギュレータGW近似A/Dコンバータカレントミラーアナログスイッチ数値微分マフィンティン半径TL431発振回路サーボPC817CUSB直流動作点解析74HC4053補間FFTBSch開発環境パラメトリック解析2ちゃんねるチョッパアンプ量子力学bzqlty電子負荷イジング模型LDA標準ロジックアセンブラ基本並進ベクトルブラべ格子単振り子熱伝導位相図TLP621キュリー温度繰り返し状態方程式MaximaVESTAスイッチト・キャパシタ相対論FETランダムウォークスピン軌道相互作用SMP六方最密充填構造抵抗不規則合金ewidthスレーターポーリング曲線GGAラプラス方程式cygwingfortranQSGW失敗談コバルト条件分岐TLP521テスタLM555Writer509TLP552格子比熱マントルデータロガー自動計測詰め回路ガイガー管ダイヤモンドQNAPMCUFXA-7020ZR過渡解析三角波UPSNE555固有値問題熱力学ブラウン運動フェルミ面awk起電力第一原理計算OpenMPfsolveubuntu最大値xcrysden最小値最適化仮想結晶近似VCA差し込みグラフスーパーセル井戸型ポテンシャル平均場近似シュレディンガー方程式FSMフラクタルOPA2277固定スピンモーメント2SC1815全エネルギー合金multiplotgnuplotc/aTeX結晶磁気異方性interp1ウィグナーザイツ胞初期値マンデルブロ集合疎行列面心立方構造fcc不純物問題非線型方程式ソルバフィルタL10構造PGA半金属二相共存SICZnOウルツ鉱構造BaO重積分クーロン散乱磁気モーメント電荷密度三次元CIF岩塩構造CapSenseノコギリ波デバイ模型ハーフメタル正規分布フォノンquantumESPRESSOルチル構造スワップ領域リジッドバンド模型edelt縮退キーボード軸ラベルグラフの分割凡例トラックボールPC不規則局所モーメント片対数グラフトランス両対数グラフCK1026MAS830L直流解析Excel円周率パラメータ・モデルヒストグラム日本語最小二乗法等価回路モデルGimp線種シンボルTS-110TS-112PIC16F785LMC662化学反応文字列specx.f入出力ifortマテリアルデザインヒストグラム確率論Realforce等高線ジバニャン方程式P-10Ubuntuナイキスト線図Crank-Nicolson法陰解法熱拡散方程式HiLAPWAACircuit連立一次方程式負帰還安定性境界条件EAGLEMBE関数フィッティング

最新コメント
リンク

にほんブログ村 その他趣味ブログ 電子工作へ