期末作業(yè)，選擇了熱力學(xué)與統(tǒng)計(jì)物理學(xué)中的內(nèi)容，一方面完成了課程考核，另一方面捎帶復(fù)習(xí)熱統(tǒng)。主要是和教材中伊辛模型?相關(guān)的熱力學(xué)與統(tǒng)計(jì)物理學(xué)中的內(nèi)容，編程畫(huà)圖分析得結(jié)論。
一，摘要
相變和臨界現(xiàn)象是一個(gè)跨學(xué)科領(lǐng)域，在連續(xù)相變的臨界點(diǎn)，熱力學(xué)函數(shù)和關(guān)聯(lián)函數(shù)呈現(xiàn)奇異性，他們來(lái)源于粒子間的相互作用。在熱統(tǒng)中，處理系統(tǒng)性問(wèn)題時(shí)往往將微觀粒子之間的相互作用作為小項(xiàng)而忽略掉，且稱這種粒子系統(tǒng)為近獨(dú)立子系，但在某些情況下，系統(tǒng)的粒子之間的相互作用不可忽略，甚至在系統(tǒng)中起著重要作用，此時(shí)就再也不能將其忽略掉了。例如物質(zhì)的鐵磁性就是粒子間較強(qiáng)相互作用引起的，處理此類系統(tǒng)問(wèn)題時(shí)便不能再忽略掉粒子間相互作用，而伊辛模型就是一種描寫(xiě)磁系統(tǒng)相變的模型，本文通過(guò)數(shù)值方法來(lái)研究伊辛模型。
二，分析與編程
對(duì)于自旋系統(tǒng)，應(yīng)用最近鄰近似法，再以一個(gè)與復(fù)雜自旋耦合效果等價(jià)的平均力場(chǎng)來(lái)代替耦合，對(duì)于粒子?i，計(jì)及其周圍共?Z?個(gè)粒子的影響，以場(chǎng)?B_e來(lái)代替這?Z?個(gè)粒子對(duì)粒子 i 的共同耦合效果。下面公式，左邊為最近鄰自旋下對(duì)粒子?i?的能量 E，右邊為以平均力場(chǎng)?B_e?作用于粒子 i?后附加的能量 E
以上引用熱統(tǒng)中的內(nèi)容來(lái)引入概念和方法，下面開(kāi)始問(wèn)題的正式解決。
對(duì)于由自旋為 1 的粒子，對(duì)于一個(gè)特定的空間取向，其自旋呈現(xiàn)只有 +1 和 -1 的變化
本課程教材中通過(guò)分析和簡(jiǎn)化，同時(shí)應(yīng)用分子場(chǎng)近似（平均力場(chǎng)）得出結(jié)論公式
<s>?為參量取某指定值時(shí)的平均自旋，J?為交換系數(shù)，k?為玻爾茲曼常數(shù)，T?為對(duì)應(yīng)溫度，z?為最近鄰粒子數(shù)目。
簡(jiǎn)單分析以上公式外加聯(lián)系以上函數(shù)的曲線特征可知，在?z?,?T?,?J?都為穩(wěn)定值時(shí)，等式兩邊的函數(shù)圖像一般存在三個(gè)交點(diǎn)，如下圖：
但是當(dāng)?tanh?函數(shù)在原點(diǎn)的斜率與直線斜率?y = <s>?相等時(shí)，便只存在一個(gè)交點(diǎn)，即等式只存在一個(gè)解 <s>=0 ，由平均場(chǎng) H'?的定義式
則 <s>=0 該情況則對(duì)應(yīng)著?H' =0 ，即磁場(chǎng)為 0，沒(méi)有自發(fā)磁化現(xiàn)象發(fā)生。由一個(gè)簡(jiǎn)單的推導(dǎo)可知，此時(shí)對(duì)應(yīng)的溫度 T = z*J/k 。
下面根據(jù)公式畫(huà)出不同溫度下曲線的交點(diǎn)情況，假設(shè) a = zJ/k ，取 a為 1（因?yàn)槲覀兊膱D像模擬只需要的是一種變化的趨勢(shì)，故常數(shù)項(xiàng) a的具體值無(wú)關(guān)緊要，取為正數(shù)即可，若要具體真實(shí)的情況，則將z，J，k代入得到a即可）
公式改造為：
當(dāng) T大于臨界值時(shí)，以此繪圖有：
當(dāng) T取為小于臨界 T值時(shí)，有關(guān)系圖：
由以上繪圖可知，當(dāng)溫度 T 大于臨界值時(shí)，直線與曲線是沒(méi)有非零交點(diǎn)的，亦即此時(shí)物質(zhì)無(wú)自發(fā)磁化現(xiàn)象發(fā)生。而當(dāng)溫度小于臨界溫度時(shí)，在一定范圍內(nèi)溫度越低，交點(diǎn)越接近?1?和? -1，即越來(lái)越接近完全的自發(fā)磁化（因?yàn)?s 最大值就是 1 ），<s>=1?即對(duì)應(yīng)著磁化方向與取定方向同向，<s>=-1?意味著與取定方向完全反向。由圖可發(fā)現(xiàn)當(dāng)溫度低到一定程度時(shí)，在此溫度以下的溫度時(shí)，都是發(fā)生完全磁化
由于原方程解析解不能求出，故用數(shù)值方法對(duì)不同溫度下的方程求解得 <s> 值
用Mathematica尋根命令：
Table[FindRoot[x == Tanh[x/a], {x, 1}], {a, 0.02, 1, 0.02}]
得方程的幾個(gè)溫度下的解：
{x -> 1.}, {x -> 1.}, {x -> 1.}, {x -> 1.}, {x -> 1.}, {x -> 1.}, {x -> 0.999999}, {x -> 0.999993}, {x -> 0.99997}, {x -> 0.999909}, {x -> 0.999774}, {x -> 0.999517}, {x -> 0.999081}, {x -> 0.998402}, {x -> 0.997414}, {x -> 0.99605}, {x -> 0.994248}, {x -> 0.991947}, {x -> 0.98909}, {x -> 0.985624}, {x -> 0.981499}, {x -> 0.976669}, {x -> 0.971089}, {x -> 0.964715}, {x -> 0.957504}, {x -> 0.949413}, {x -> 0.940398}, {x -> 0.930412}, {x -> 0.919408}, {x -> 0.907332}, {x -> 0.894128}, {x -> 0.879732}, {x -> 0.864073}, {x -> 0.847072}, {x -> 0.828635}, {x -> 0.808656}, {x -> 0.78701}, {x -> 0.763548}, {x -> 0.738089}, {x -> 0.710412}, {x -> 0.680241}, {x -> 0.64722}, {x -> 0.610884}, {x -> 0.570591}, {x -> 0.52543}, {x -> 0.474002}, {x -> 0.413976}, {x -> 0.340829}, {x -> 0.242983}, {x -> 2.25808*10^-8}
畫(huà)圖有：
由上圖可見(jiàn)，在一定范圍溫度內(nèi)，<s> 保持完全磁化狀態(tài)，當(dāng)溫度超過(guò)某值后開(kāi)始下降，且下降趨勢(shì)越來(lái)越快，當(dāng)達(dá)到臨界溫度（對(duì)應(yīng)橫坐標(biāo)為 1.0?）后自發(fā)磁化現(xiàn)象消失。
用蒙特卡洛方法對(duì)磁系統(tǒng)的研究
本部分由教材內(nèi)容直接應(yīng)用，簡(jiǎn)略描述為：自旋 s_i?受其周圍諸最近鄰自旋?s_j?的影響，會(huì)發(fā)生取向反轉(zhuǎn)，對(duì)于二維方格子情形，因?yàn)樽孕谥辉诮o定方向上取值 ±1?，所以之前的能量公式由矢量式退化為標(biāo)量式：
選取一包含方陣結(jié)構(gòu)一定數(shù)量自旋粒子的方格子，粒子 i 受周圍最近鄰作用：
應(yīng)用蒙特卡洛方法的原理??是，先給定一個(gè)格子內(nèi)部各自旋的初始分布，以粒子 i 為例，當(dāng)其自旋反轉(zhuǎn)時(shí)，則顯然能量 E 的值會(huì)變化，自旋反轉(zhuǎn)前后能量的變化值：
對(duì)上式的解釋：E'?即反轉(zhuǎn)后能量，E_0 則為原能量。將?ΔE?代入到玻爾茲曼因子?e^(-ΔE/kT)中。改寫(xiě)為
由隨機(jī)數(shù)發(fā)生器產(chǎn)生 [0,1]?之間的隨機(jī)數(shù)。ΔE > 0?時(shí)，當(dāng)隨機(jī)數(shù)小于上式值時(shí) S_?i?反轉(zhuǎn)，當(dāng)隨機(jī)數(shù)大于上式值時(shí)?S_?i?不反轉(zhuǎn)；而當(dāng)?ΔE ≤?0 時(shí)，S_?i反轉(zhuǎn)，邊界處應(yīng)用周期性邊界條件，以此規(guī)則編程即可。編程時(shí)以 t =?2J/kT 為自變量。
t = 5; l = Length[s]; list = {};
Do[s = Table[{1, -1, 1, 1, -1, 1, -1, -1, 1, -1, 1, -1, 1, -1, -1, 1,1, -1, 1, -1, 1, -1, 1, -1, 1, -1, 1, 1, -1, 1}, 30];
? sum = Total[Flatten[s]] - 100;
? While[Abs[Total[Flatten[s]] - sum] >= 1, sum = Total[Flatten[s]];?
For[i = 1, i <= 30, i++,
? ? For[j = 1, j <= 30, j++,
? ? If[i == 1 && j >= 2 && j <= l - 1,
? ? ? dE = 2 s[[i,j]]*(s[[l, j]] + s[[i + 1, j]] + s[[i, j - 1]] + s[[i, j + 1]]), Null];
? ? If[i == l && j >= 2 && j <= l - 1, dE = 2 s[[i, j]]*(s[[i - 1, j]] + s[[1, j]] + s[[i, j - 1]] +s[[i, j + 1]]), Null];
? ? If[j == 1 && i >= 2 && i <= l - 1,dE = 2 s[[i, j]]*(s[[i - 1, j]] + s[[i + 1, j]] + s[[i, l]] + s[[i, j + 1]]), Null];
? ? If[j == l && i >= 2 && i <= l - 1,dE = 2 s[[i, j]]*(s[[i - 1, j]] + s[[i + 1, j]] + s[[i, j - 1]] +s[[i, 1]]), Null];
? ? If[i == 1 && j == 1,dE = 2 s[[i, j]]*(s[[l, j]] + s[[i + 1, j]] + s[[i, l]] + s[[i, j + 1]]), Null];
? ? If[i == 1 && j == l,dE = 2 s[[i, j]]*(s[[l, j]] + s[[i + 1, j]] + s[[i, j - 1]] + s[[i, 1]]),Null];
? ? If[i == l && j == l, dE = 2 s[[i, j]]*(s[[i - 1, j]] + s[[1, j]] + s[[i, j - 1]] + s[[i, 1]]), Null];
? ? If[i == l && j == 1, dE = 2 s[[i, j]]*(s[[i - 1, j]] + s[[1, j]] + s[[i, l]] + s[[i, j + 1]]), Null];
? ? If[i > 1 && j > 1 && i < l && j < l, dE = 2 s[[i, j]]*(s[[i - 1, j]] + s[[i + 1, j]] + s[[i, j - 1]] +
? ? ? ? ? s[[i, j + 1]]), Null];
? ? If[dE <= 0 || E^(-t*dE) >= RandomReal[], s[[i, j]] = -s[[i, j]],? Null]
? ? ]]]; AppendTo[list, sum/l^2], 100];
ListPlot[list]
以上編程中取了一個(gè)30*30的自旋方格子，運(yùn)行結(jié)果：
t=0.1
t=0.5
t=1
t=2
可以看到隨著t的增大（因?yàn)?i> t 與溫度 T 是反比關(guān)系，所以隨著 t 的增大亦即隨著溫度 T 的減?。?，自發(fā)磁化現(xiàn)象越來(lái)越明顯。在溫度T較高時(shí)震蕩較大，且是圍繞著無(wú)自發(fā)磁化震蕩，而隨著溫度T的減小，自發(fā)磁化開(kāi)始出現(xiàn)，且較明顯。由上圖即可直觀地發(fā)現(xiàn)溫度與自發(fā)磁化（方格子自選統(tǒng)計(jì)量）的關(guān)系。以上即運(yùn)用了蒙特卡洛方法研究伊辛模型的模擬結(jié)果。
三，總結(jié)
本作業(yè)主要做了兩件事，一個(gè)是應(yīng)用一般數(shù)值方法?研究磁系統(tǒng)的自發(fā)磁化強(qiáng)度隨溫度的變化，二是應(yīng)用所謂的蒙特卡洛方法?研究自旋方格子內(nèi)的磁化強(qiáng)度隨溫度的變化
通過(guò)模擬得到了一些結(jié)論，因?yàn)橐列聊Ｐ偷木_解難以求解，而計(jì)算機(jī)的告訴數(shù)值運(yùn)算能力便應(yīng)用解決了這一瓶頸。從以上模擬中，能看出自旋間相互作用隨外界溫度的變化情況，符合理論上的定性描述，而一些關(guān)鍵相變電也被數(shù)值方法描述了出來(lái)，總而言之，此次數(shù)值模擬是驗(yàn)證了理論的正確性，也能看出，數(shù)值模擬能彌補(bǔ)解析方法的不足，帶來(lái)一些直觀的圖像。
四，說(shuō)明
參考文獻(xiàn)有：
本課程教材
武漢大學(xué)出版社? 胡承正《熱力學(xué)與統(tǒng)計(jì)物理學(xué)》
二維Ising模型臨界相變的Monte-Carlo數(shù)值模擬（《安徽大學(xué)學(xué)報(bào)：自然科學(xué)版》2008年 第3期 | 李曉寒 王宗笠 寧旭 ）
作業(yè)最后一次更新于2017/12/28，并于該日將pdf文件郵件發(fā)至老師郵箱hcai@outlook.com
