MATLABで簡単なシミュレーションをしよう

はじめに

大学，大学院では数値計算ソフトとしてMATLABを使う方も多いのではないでしょうか．
やはり最近は，python人気ということもありMATLABが見向かれなくなっているように感じますが，数値計算全般ではコードが非常に組みやすく，pythonに比べても構文が簡単に書けることも多いと思っています．
今回は，MATLABを研究で使い始める人や初めて使う人向けに，具体的な使用例を紹介したいと思います．

量子力学

MATLABの使用例として，「線形代数の計算」と「結果の可視化」という最も基礎的な例を取り上げます．その対象として，系のエネルギーをハミルトニアンで記述できる量子力学を例題として挙げようと思います．

ハミルトニアン

量子力学において，系のエネルギーを記述するのはハミルトニアンでした．ハミルトニアンの固有値がエネルギーです．いま

「電子スピンに外部磁場を印加した際に，その磁場強度と磁場角度によってどのようなエネルギーをとるのか？」

を可視化します．ここではMATLABの練習が主題なので物理的意味は興味なければ追求する必要はありません．
（※スピンの量子化軸を$z$軸方向にとり，$z$方向とのなす角度を$\theta$と定義します）

ハミルトニアンは

H = DS_z^2\,+\,\gamma B_0S_{B_0} \\
S_{B_0} = S_z\cos\theta + S_x\sin\theta

で表します．問題は，磁場角度θと磁場Bです．これらを振って，Hの固有値を求めましょう（他のパラメータは全て定数です）．

tripletstate.m

clc
clear

Sx = 1/sqrt(2)*[0,1,0;1,0,1;0,1,0];
Sz = [1,0,0;0,0,0;0,0,-1];
Ns = size(Sz,1); 

% physical quantity
gamma = 2.8e+6;% Hz/G
Dgs = 2.87e+9;% GHz
Hdgs = Dgs*Sz.^2;

B0 = linspace(0,1200,1e+2); %等差数列を作る
Nb0 = length(B0);

theta = linspace(0,pi/2,1e+2);
Ntheta = length(theta);

Energy = zeros(Nb0,Ntheta,Ns); %３次元のゼロの配列
for refth = 1:Ntheta
    Sb0 = Sz*cos(theta(refth))+Sx*sin(theta(refth));
    for refb0 = 1:Nb0
        Hzeeman = gamma*B0(refb0)*Sb0;

        Htot = Hdgs+Hzeeman;

        E = eig(Htot);
        Energy(refb0,refth,1) = E(1);
        Energy(refb0,refth,2) = E(2);
        Energy(refb0,refth,3) = E(3);
    end
end

%グラフの描画
figure  
hold on
box on
for refl = 1:Ns
    surf(B0, theta/pi*180, Energy(:,:,refl)'*1e-9);
end
hold off

xlabel('B_0 (G)');
ylabel('\theta (deg)');
zlabel('Energy (GHz)');
shading('flat');
colorbar;

コードはこれだけです．
始めに，物理定数$D$と行列$S_x, S_z$を定義します．ここではパウリ行列を用いています．
for文直前には，計算領域の確保を行います．必要なのは，2変数に対する結果なので計3次元の領域が必要になり，zerosを用いて0行列を作製しておく必要があります（ここに順次，計算結果を格納します）．

あとはfor文中で，$\theta$と$B_0$の各要素に対するハミルトニアンを計算します（要素数を予め，Nb0, Nthetaとして取得しておくのがよいです）．

ハミルトニアンを計算し，最後にeigを用いて固有値が算出されます！便利．

当然，各々の$\theta$と$B_0$に対して計算するのでfor文内でEnergy(refb0,refth,1)などに順次入れていく必要がありますね．

グラフの描画

計算結果を踏まえて，グラフを描画します．
3次元にはsurfが便利です．また，shading('flat')をいれることで滑らかな曲線を描けます．

結果

このハミルトニアンを解いた結果，以下のようなグラフが得られました．

大学を卒業後に，MATLABは使えなくなってしまった（笑）ので，今回はOctaveで出力しています．しかしコードは全く同じで大丈夫です．3D図も描けます．
物理的な意味は詳述しませんが，ここでは磁場と角度を変数としてEnergy（縦軸）の準位が計算できました．固有値をしっかり表現できたようです．

おわりに

今回はハミルトニアンの固有値を得て，それを可視化しました．
その過程では，

１．行列計算
２．固有値計算
３．グラフの出力

を行う大変基本的なものですが，計算結果の格納，for文という非常に重要な要素も使われています．

私はコードを書く勉強を始めたとき，何からやればいいかわからず最初は先輩に書いてもらいました．最初は真似から入るのが良いと思います．
とくに，物理を始めた学生方は大学の計算環境を用いていかがでしょうか．