Scilabで乱数の生成

Scilabで楽しむ確率論(PDF)に倣って乱数を使ったモンテカルロ法を勉強します。


004_2013092414520065a.png
Fig.1: 一様分布のばらつきとヒストグラム。



これまでたくさんの常微分方程式タグの付いたエントリで微分方程式による物理現象のモデル化(PDF)の計算をScilabで行ってきました。
次の章からは統計熱力学を扱うことになるのですが、ねがてぃぶろぐではその前にScilabでの乱数の使い方を勉強するためにScilabで楽しむ確率論(PDF)に行こうと思います。

一様分布の乱数


シミュレーションや数値計算に乱数を使う方法を総じてモンテカルロ法と呼びます。回路シミュレータLTspiceにも乱数を使ったモンテカルロ解析が実装されているのでした(参考:LTspiceでモンテカルロ解析)。
LTspiceモンテカルロ解析の定数分布 その1ではLTspiceのMC関数で使われる乱数の分布を調べました。

Scilabで乱数を使うにはrandを使います。

plot(rand(1,10000));


とすると、乱数がプロットされます。


001_201309241451268ad.png
Fig.2: 横軸が繰り返し番号で縦軸が発生された乱数の値。


グラフより0から1までの値が無秩序に並んでいることが分かります。

次にhistplotを使ってヒストグラムを書きます。LTspiceのときは回りくどい方法を使った記憶がありますが、Scilabでは簡単にプロットできます。

histplot(100,rand(1,10000));



002_20130924145201cc4.png
Fig.3: 一様乱数のヒストグラム


乱数が一様分布であることが一目で分かります。

三角分布の乱数


一様分布を加減算すると三角分布になるのでした。(参考:LTspiceモンテカルロ解析の定数分布 その6)

Scilabでも同じ計算をやってみます。

clear;

// 三角分布の発生
X = 0.5 * (rand(1,10000) + rand(1,10000));

// グラフのプロット
subplot(2,1,1);
plot(X);
// ヒストグラム
subplot(2,1,2);
histplot(100,X);



003_201309241452005d3.png
Fig.4: 三角分布


確かに三角分布が得られました。

関連エントリ




参考URL




付録


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


参考文献/使用機器




フィードバック



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

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


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


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

tag: Scilab モンテカルロ解析 乱数 ヒストグラム確率論 

AkaiKKRで鉄のキュリー温度

強磁性体のキュリー温度は、平均場近似から以下のように求めることが出来ます。

T_c = \frac{2}{3ck_B}\Delta E

(Tc: キュリー温度, c: 磁性原子の濃度 (0 < c ≦ 1), kB: ボルツマン定数, ΔE: 常磁性状態と強磁性状態の間のエネルギー差)

計算機マテリアルデザイン入門 (大阪大学新世紀レクチャー)密度汎関数法の発展 -マテリアルデザインへの応用では、平均場近似とAkaiKKR(Machikaneyama)による全エネルギー計算による強磁性体のキュリー温度の見積もりが紹介されています。

今回は、これらに倣って鉄のキュリー温度を計算しました。計算結果は Tc = 1260 (K) となり、鉄のキュリー温度の実測値である Tc = 1043 (K) と近い値が得られました。

001_20130924035510a8c.png 002_201309240355099cc.png


AkaiKKRと平均場近似によるキュリー温度の計算


計算機マテリアルデザイン入門 (大阪大学新世紀レクチャー)によると、強磁性体のキュリー温度は平均場近似を用いて以下のように求めることが出来ます。

T_c = \frac{2}{3ck_B}\Delta E

(Tc: キュリー温度, c: 磁性原子の濃度 (0 < c ≦ 1), kB: ボルツマン定数, ΔE: 常磁性状態と強磁性状態の間のエネルギー差)

強磁性状態のエネルギーは、AkaiKKRでニッケル・鉄・コバルトで行った様に簡単に計算できます。一方で、常磁性状態の計算は非磁性(nmag)の計算とは異なります。

(なのでAkaiKKRで鉄の安定相と格子定数の非磁性の計算をδ鉄と呼ぶのは間違いだったということです。すみません。)

常磁性状態は、それぞれの原子の置ける磁気モーメントがランダムな向きを向いている状態です。しかしながら、第一原理計算でこの状態を再現するのは難しいようです。

そこでAkaiKKRでは、常磁性状態の代わりに不規則局所モーメント状態の計算を行います。これは上向きの磁気モーメントをもつ原子と下向きの磁気モーメントを持つ原子が半分ずつ不規則合金となったもので、局所的には磁気モーメントを持ちながら、金属結晶全体ではモーメントが無いという特徴を持っています。

今回は『平均場近似』と『常磁性状態の変わりに不規則局所モーメント状態を計算する』という2つの近似の下に鉄のキュリー温度を計算します。

fmg.fのコンパイル


CygwinでAkaiKKR(Machikaneyama)ではAkaiKKR本体をcygwinのg77でコンパイルしました。AkaiKKRには不規則局所モーメント状態のための初期ポテンシャルデータを作成するための補助プログラムとしてfmg.fが付属していますので、これをコンパイルします。

計算機マテリアルデザイン入門 (大阪大学新世紀レクチャー)のP255に書いてある通りCygwinの端末上でutilフォルダに移動した後
~/cpa2002v009c/util> f77 -o fmg fmg.f

とタイプします。コンパイルが成功すればfmg.exeという実行ファイルが出来ているはずです。

計算手順


第一原理計算の部分をまとめて実行するために、あらかじめ必要なファイルを全て用意しておきましょう。
準備する入力ファイルは下記の3つです。
  • 強磁性状態のための入力ファイル: fefmg.in
  • 不規則局所モーメント状態のための入力ファイル: felmd.in
  • fmg.exeのための入力ファイル: fefmg


キュリー温度を求めるためには、2つの状態の全エネルギーを求めるだけでよいのですが、今回は状態密度も同時に計算できる入力ファイルを準備しました。
各入力ファイルの準備が完了したら、いよいよAkaiKKRを用いた第一原理計算を実行します。
手順は以下のようになります。

  1. 強磁性状態の計算
  2. fmg.exeを利用した強磁性状態のポテンシャルファイルからの不規則局所モーメント状態のための初期ポテンシャル作成
  3. 不規則局所モーメント状態の計算


上記の手順を人間が行っても良いのですが、一気にやってくれるシェルスクリプトFe.shを用意しました。

Curie/─┬─in/─┬─fefmg.in
│ └─felmd.in
├─out/
├─data/
├─Fe.sh
└─fefmg


上記のようなディレクトリ構成としてFe.shをCygwin端末上で実行します。

キュリー温度


3度目の掲示になりますが、下記が平均場近似によるキュリー温度の算出式です。

T_c = \frac{2}{3ck_B}\Delta E

(Tc: キュリー温度, c: 磁性原子の濃度 (0 < c ≦ 1), kB: ボルツマン定数, ΔE: 常磁性状態と強磁性状態の間のエネルギー差)

ここで磁性原子の濃度はc=1です。
AkaiKKRによる全エネルギーはRyの単位で出力されるのでJへ変換する必要があります。

1 (Ry) = 2.179 872×10-18 (J)

またボルツマン定数は

kB = 1.3806488×10-23 (J/K)

です。

Efmg = -2522.8176206 (Ry)
Elmd = -2522.8055978 (Ry)

なので、求められた全エネルギーからキュリー温度は Tc = 1259.194 (K) と計算されました。この値は鉄のキュリー温度の実測値である Tc = 1043 (K) と近い値です。

状態密度


ついでに計算を行った強磁性鉄と不規則局所モーメント鉄の状態密度を示します。


001_20130924035510a8c.png

Fig.1: 常磁性状態の状態密度

002_201309240355099cc.png

Fig.2: 不規則局所モーメント状態の状態密度


状態密度のプロットは以下のファイルによって行いました。



Appendix: 改行コードの問題


WindowsでCygwinを使う上での特有の事ですが、改行コードの問題があります。
Windowsでは普通改行コードにCR+LFを用います。ところがLinuxではLFを利用します。

今回利用したシェルスクリプトFe.shとfmg.exeのための入力ファイルfefmgはどちらもLinux流のLFの改行コードで無いと正常に動作しません。

エラーが出た場合は確認してみてください。

関連エントリ




参考URL




付録


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



参考文献/使用機器





フィードバック



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

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


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


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

tag: AkaiKKR machikaneyama KKR CPA 強磁性 キュリー温度 平均場近似 不規則局所モーメント 状態密度 DOS 

Scilabのコンソール画面とデータのやり取り

計算物理学 基礎編の最初のプログラムである円の面積の計算をScilabで作成し、Scilabのプログラムとコンソール画面とのデータのやり取りについて練習します。(参考:area.f, area.c)


円の半径を入力させて、それに応じた円の面積を計算するScilabプログラムはaria_sce.txt
です。

clear;

// 半径rの入力を促す
r = input("Enter the radius of a circle: ");
// 円の面積を計算
a = %pi * r ^ 2;

// 円の半径と面積を表示する
disp("r = " + string(r) + ", " + "A = " + string(a));


数値の入力を促す命令がinputです。

文字列や数値をコンソール画面に出力する命令がdispです。

disp("Hello, Scilab!");


変数を文字列に変換する命令はstringです。

文字列をつなげるのは、数値を足し算するのと同様に + 記号を使います。普通の足し算と同じ記号を使うのは紛らわしいなと感じる場合はstrcatを使うのが良いと思います。

参考URL




付録


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


参考文献/使用機器




フィードバック



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

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


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


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

tag: Scilab 文字列 入出力 

Scilabでガウス型波束の散乱

Scilabを用いて初期波動関数がガウス型をしている波束が、ポテンシャルの谷によって散乱される様子をシミュレーションしました。

001.gif

Fig.1: ガウス型波束がポテンシャルの谷間によって散乱される様子を表したアニメーション。赤が波動関数の絶対値、緑が波動関数の実部を表す。青はポテンシャルの形状。



偏微分方程式


ねがてぃぶろぐでは微分方程式による物理現象のモデル化(PDF)で紹介されている常微分方程式をScilabで計算してきました。
次のテーマである波束の運動の問題は、時間微分だけの常微分方程式ではなく、時間と空間の両方の微分を含む偏微分方程式です。

ScilabやOctaveには偏微分方程式を解くための特別な関数は用意されていません。そこで、偏微分方程式を解くためには微分を差分に置き換えた計算が必要になります。Scilabで一次元井戸型ポテンシャルScilabで規則ポテンシャルと縮退では位置の微分であるエネルギー演算子を含むハミルトニアンを行列で表すことによって、固有値問題を解くための命令であるspecを利用することにより計算しました。

今回はScilabで一次元井戸型ポテンシャルと同様にハミルトニアンを行列であらわし、常微分方程式ソルバodeと組み合わせることにより、時間に依存するシュレディンガー方程式を解いてガウス型波束が散乱される様子をシミュレーションします。


ガウス型波束の運動


シッフの量子力学の本(上巻 下巻)によるとx軸の正の方向に運動量を持つ波束の波動関数は以下のように表されます。

\psi(x) = \frac{1}{\sqrt[4]{2 \pi (\Delta x)^2}} \exp\left(- \frac{(x-\langle x \rangle)^2}{4(\Delta x)^2}+\frac{i\langle p \rangle}{\hbar}x \right)

今回は、このガウス型波束を初期値に持つ波動関数がポテンシャルの谷に散乱される様子をScilabを用いてシミュレーションします。


時間に依存するシュレディンガー方程式


時間に依存するシュレディンガー方程式は、以下の様に表されます。

i \hbar \frac{\partial \psi (x,t)}{\partial t} = \left(- \frac{\hbar}{2m}\frac{\partial^2}{\partial x^2} + V(x) \right)\psi(x,t)

時間に依存しないシュレディンガー方程式の場合と同様に

\hbar = 1, m = \frac{1}{2}

とおくと、ハミルトニアンHは以下のようにかけます。

H = - \frac{\partial^2}{\partial x^2} + V(x)

この段階で時間に依存するシュレディンガー方程式は以下のように書くことが出来ます。

i\frac{\partial \psi}{\partial t} = H \psi

ここまでくればScilabの常微分方程式ソルバodeで解けることが分かります。
dψ/dt = ...の形にすると

\frac{\partial \psi}{\partial t} = -i H \psi

となり、これが今回の解くべき方程式となります。

微分方程式による物理現象のモデル化(PDF)によるとOctaveの常微分方程式ソルバであるlsodeは複素数に対応していないとのことですが、Scilabのodeは複素数に対応しているので、今回はodeを利用してシュレディンガー方程式を解きます。


疎行列


Scilabで一次元井戸型ポテンシャルで書いたとおり、ハミルトニアンを行列として書き下すと、以下のようにほとんどの成分がゼロとなります。

\begin{eqnarray*} H &=& - \frac{\mathrm{d}^2}{\mathrm{d}x^2} + V \\ &=& \left( \begin{array}{ccccccc} \frac{2}{h^2} + v_1 & -\frac{1}{h^2} & 0 & \hdots & & \hdots & 0 \\ -\frac{1}{h^2} & \frac{2}{h^2} + v_2 & -\frac{1}{h^2} & 0 & \hdots & & \vdots \\ 0 & -\frac{1}{h^2} & \frac{2}{h^2} + v_3 & -\frac{1}{h^2} & 0 & \hdots & \\ \vdots & \ddots & \ddots & \ddots & \ddots & \ddots & \vdots \\ & \hdots & 0 & -\frac{1}{h^2} & \frac{2}{h^2} + v_{n-2} & -\frac{1}{h^2} & 0 \\ \vdots &  & \hdots & 0 & -\frac{1}{h^2} & \frac{2}{h^2} + v_{n-1} & -\frac{1}{h^2} \\ 0 & \hdots &  & \hdots & 0 & -\frac{1}{h^2} & \frac{2}{h^2} + v_n \end{array} \right) \end{eqnarray*}

このようにゼロでない要素が全体の中でごくわずかしかない行列を疎行列と呼びます。Scilabでは疎行列を扱う特別なデータ型があり、普通に定義した行列に対して

sp = sparse(X);


とすることによって疎行列へ変換することが出来ます。(参考:sparse)

今回の計算のハミルトニアンのようにゼロとなる要素が多い場合は、疎行列型へ変換することにより計算速度が劇的に速くなります。


数値計算


微分方程式による物理現象のモデル化(PDF)のように100stepごとに700ステップまで表示するScilabのプログラムがschiff-gv_sce.txtです。

clear;

vw = 0.064; // ポテンシャルの幅
vh = - (70.7 * %pi) ^ 2; // ポテンシャルの深さ
dim = 500; // 空間の分割数
h = 2.0 / dim; // 空間の刻み幅

// *** 位置ベクトル ***
X = linspace(-1,1,dim + 1);

// *** ハミルトニアンの定義 ***
// エネルギー演算子
vd = 2 / h ^ 2;
vdt = - 1 / h ^ 2;
SD = diag(vdt * ones(1,dim),1);
D = diag(vd * ones(1,dim + 1));
H = SD + D + SD.';
// ポテンシャル
V = zeros(H);
for i = 1:dim + 1 do
if abs(X(i)) < vw / 2 then
V(i,i) = vh;
end
end
// 全ハミルトニアン
H = sparse(H + V); // 疎行列へ変換するほうが実行速度が速い

// *** 常微分方程式の定義 ***
function dphi = schiff(t,phi)
dphi = - %i * H * phi;
endfunction

// *** 波束の初期設定 ***
xe = - 0.3; // 初期波束の位置の期待値(波束の中心位置)
dx = 0.035; // 初期波束の空間的広がり
pe = 50 * %pi; // 初期波束の運動量の期待値
phi0 = zeros(dim + 1,1);
for k = 1:dim + 1 do
// 初期波束
phi0(k) = exp(-(X(k) - xe) ^ 2 / (4 * dx ^ 2) + %i * pe * X(k)) ..
/ (2 * %pi * (dx) ^ 2) ^ 0.25;
end

// *** 時間刻みの設定 ***
dt = h ^ 2 / 4;
nstep = 700;
itvl = 100;
tf = dt * itvl;
T = [0:dt:tf];

// *** プロット ***
// ポテンシャルのプロット
subplot(3,3,1);
plot(X,diag(V));
// 目盛りを非表示にする
g = gca();
g.axes_visible = 'off';
zoom_rect([-0.5,-5,0.5,5]);
// 初期値のプロット
subplot(3,3,2);
plot(X,diag(V)); // ポテンシャルのプロット
plot(X,real(phi0),'-g'); // 波動関数の実部のプロット
plot(X,abs(phi0),'-r'); // 波動関数の絶対値のプロット
// ステップ数を表示
xstring(0.1,3,"step = 0");
// 目盛りを非表示にする
g = gca();
g.axes_visible = 'off';
zoom_rect([-0.5,-5,0.5,5]);

for k = 1:nstep/itvl do
// 微分方程式の数値解
phi = ode(phi0,0,T,schiff);
// ポテンシャルと波動関数のプロット
subplot(3,3,k + 2);
plot(X,diag(V)); // ポテンシャルのプロット
plot(X,real(phi(:,itvl)),'-g'); // 波動関数の実部のプロット
plot(X,abs(phi(:,itvl)),'-r'); // 波動関数の絶対値のプロット
// ステップ数を表示
stepnum = string(k * itvl);
xstring(0.1,3,strcat(["step = ",stepnum]));
// 目盛りを非表示にする
g = gca();
g.axes_visible = 'off';
zoom_rect([-0.5,-5,0.5,5]);
// 最後の波動関数を次の初期値に設定
phi0 = phi(:,itvl);
end


Fig.1のGIFアニメーションはschiff-gv-gif_sce.txtで作成したGIF画像をGimpやImagemagickで結合したものです。
Scilabの計算結果をGIFアニメーションにする方法はのちのち別のエントリにまとめようと思っています。

GIF動画に関連する情報は以下のリストにあります。


関連エントリ




参考URL




付録


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


参考文献/使用機器




フィードバック



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

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


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


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

tag: Scilab ode 常微分方程式 偏微分方程式 疎行列 量子力学 

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

LTspiceAkaiKKRScilabmachikaneyamaKKRPSoCOPアンプPICCPA強磁性モンテカルロ解析常微分方程式トランジスタode状態密度インターフェースDOSPDS5022ecaljスイッチング回路定電流半導体シェルスクリプトレベルシフト乱数HP6632A温度解析ブレッドボードI2CR6452A分散関係トランジスタ技術可変抵抗確率論数値積分反強磁性セミナー非線形方程式ソルバ絶縁バンドギャップ熱設計偏微分方程式バンド構造GW近似カオス三端子レギュレータLEDフォトカプラシュミットトリガISO-I2CA/DコンバータLM358USBカレントミラーTL431マフィンティン半径PC817C数値微分アナログスイッチ発振回路サーボ直流動作点解析74HC40532ちゃんねる標準ロジックチョッパアンプLDAアセンブラFFTbzqltyイジング模型ブラべ格子開発環境補間量子力学電子負荷BSchパラメトリック解析単振り子基本並進ベクトル熱伝導繰り返しGGAMaximaTLP621ewidthSMP相対論抵抗位相図ランダムウォークスピン軌道相互作用六方最密充填構造不規則合金FETコバルト失敗談QSGWcygwinスレーターポーリング曲線スイッチト・キャパシタラプラス方程式gfortranキュリー温度状態方程式条件分岐格子比熱TLP552LM555TLP521三角波NE555過渡解析FXA-7020ZRWriter509テスタ詰め回路MCUマントルダイヤモンドQNAPデータロガーガイガー管自動計測UPS井戸型ポテンシャルawk第一原理計算仮想結晶近似ブラウン運動差し込みグラフ平均場近似fsolve起電力熱力学OpenMPスーパーセル固有値問題最適化最小値VCAシュレディンガー方程式VESTAubuntu最大値面心立方構造PGAOPA2277L10構造非線型方程式ソルバ2SC1815fccフェルミ面等高線ジバニャン方程式ヒストグラム確率論マテリアルデザイン正規分布結晶磁気異方性interp1フィルタ初期値ウィグナーザイツ胞c/aルチル構造岩塩構造スワップ領域リジッドバンド模型edeltBaOウルツ鉱構造重積分SIC二相共存ZnOquantumESPRESSOCapSensegnuplotmultiplot全エネルギー固定スピンモーメントFSM合金ノコギリ波フォノンデバイ模型ハーフメタル半金属TeXifortTS-110不規則局所モーメントTS-112等価回路モデルパラメータ・モデルヒストグラムExcel円周率GimpトラックボールPC直流解析入出力文字列マンデルブロ集合キーボードフラクタル化学反応三次元Realforce縮退日本語最小二乗法関数フィッティング疎行列シンボル線種ナイキスト線図陰解法負帰還安定性熱拡散方程式EAGLECrank-Nicolson法連立一次方程式P-10クーロン散乱Ubuntu境界条件MBEHiLAPW軸ラベルトランスCK1026MAS830L凡例PIC16F785LMC662AACircuit両対数グラフ片対数グラフグラフの分割specx.f

最新コメント
リンク

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