数值代数上机作业第七次
二维Ising模型的相变现象
谭昌汇00401146
北京大学数学科学学院科学与工程计算系，100871
完成于2007年6月19日
摘要(Abstract)
本报告主要讨论了二维Ising模型中的内能和比热容随温度变化的相变现
象。通过对内能和比热容图像的突变情况，可以确定相变的临界温度 c。
为解决二维Ising模型，本文使用了Monte-Carlo方法，利用Metropolis算
法生成一组随机向量，从而对二维Ising模型进行求解。
关键词：二维Ising模型，相变现象，临界温度，Monte-Carlo方法
1 介绍(Introduction)
在当代的物理学中，相变问题是一个很前沿的课题，它在很多实际的问题
中都起到了重要的作用。而Ising模型是描述铁磁相变的最具代表性的模型。
对于二维的Ising模型，在理论分析上有重要的成果，即可以求得严格解。
但尽管如此，并不能给出物理图像。因此，使用计算机处理Ising模型成为了
一种十分重要的手段。Monte-Carlo方法就是其中的一种。它可以求出模型的
解，并且还可以给出具体的图像与涨落。下面的所有讨论都是围绕使用Monte-
Carlo方法求解二维Ising模型的相变问题而展开的。
2 相变问题与Ising模型
(The problem of phase change and Ising model)
2.1 相变问题(The problem of phase change)
这里所谓的相变，就是指内能、比热容等参数随温度变化产生的突变。这
些参数的突变反映了铁磁系统丧失磁性的相变过程。解决相变问题的关键就是
将不同温度下的内能和比热容的图像展示出来，从而找出临界温度。
2.2 Ising模型初探(The first approach to Ising model)
在统计物理中，人们用晶格模型研究磁性材料的相变，每个格点代表一个
电子，它有两种自旋的方式，分别记为1和-1。那么，整个系统的内能和比热容
可以有以下公式导出。
UM =
1
M
X
expf H( )g
ZM
H( ) ,
1
M
hH( )i
1
CM =
2
M
f
X
H2( )
expf H( )g
ZM
[
X
H( )
expf H( )g
ZM
]g =
2
M
[hH2( )ihH( )i2]
其中，H( ) = J
P
i j是能量函数，满足Gibbs分布。具体来说，就是相邻项
的乘积之和。比如对于下图，S1状态只需要和与其相邻的S2至S5状态乘积加和
即可。
S2 S1
S5
S3
S4
2.3 遍历性(Ergodicity)
在前面的公式中提到了hH( )i是能量状态的空间平均，是比较难于计算
的。然而，我们可以使用Monte-Carlo方法的思想，进行下面的估计。
hH( )i
1
N
XN
i=1
H( (i)):
其中，右边的序列f (i)gN
i=1是马氏链。左式为空间平均，右式为时间平均，这种
空间平均等于时间平均的结果称之为遍历性。
有遍历性的保证，即可用时间序列估计空间平均，而由Monte-Carlo方法思
想，只要适当选取f (i)gN
i=1，就可以对时间平均进行估计，最终得到上面的公式
结论。
2.4 Metropolis算法(Metropolis algorithm)
剩下的工作就是去生成一组合适的f (i)gN
i=1，我们可以使用Metropolis算法来
达到这个目的。这里仅给出算法，不进行进一步的说明，更多的信息可参见参
考文献[1]。
[算法] Metropolis算法
步1：根据相应的预选规则由 (n)产生预选态 0.
步2：计算A = minf1; exp( H)g，其中 H = H( 0) H( ).
步3：生成均匀分布的随机数r v U [0; 1].
步4：如果r A，则 n+1 = 0；否则 n+1 = n，转步1.
所谓的预选态，一般有两种，一是等可能预选，即在非 (n)状态中随机
选择一个状态。另一种是单点变化预选，变化 n的一个格点取值，总共M（格
点数）种不同情况，随机取一个即可。
使用Metropolis算法可直接调用metropolis.m文件即可。其中Hsigma参数指
的是H( )，是为了加快运算速度而加入的（通过后面的分析可以发现这样的优
化很有必要），读入时直接写作f uncH(sigma)即可。
2
2.5 小结(Brief summary)
总结一下，讨论二维Ising模型的相变现象，首先使用Metropolis算法生成时
齐马氏链，然后利用遍历性和Monte-Carlo算法估计空间平均，最后根据公式求
得不同温度下的内能和比热容，作出图形，找到突跃点，即相变的临界点。
理论已经完备，下面的工作就是进行实际的数值实验。
3 数值实验(Numerical experiment)
3.1 生成能量函数(Creating power function)
前面已经大致提到了能量函数选取Gibbs分布，即紧邻格点的求和。对于边
界条件，需要做一些必要的讨论。
我们选取的边界条件是周期条件，即对M M格点来说，第一列（行）和
最后一列（行）的自旋方向是对应相等的。要生成这样的分布，只需要先生成
一个任意的(M 1) (M 1)的分布，再扩张得到一个满足条件的Gibbs分布。具
体可见下表。
0BBBBBBBBBB@
11 1(M1)
:::
: : :
:::
(M1)1 (M1)(M1)
1CCCCCCCCCCA
=)
0BBBBBBBBBBBBBBB@
11 1(M1) 11
:::
: : :
:::
:::
(M1)1 (M1)(M1) (M1)1
11 1(M1) 11
1CCCCCCCCCCCCCCCA
对(M 1) (M 1)的分布使用Metropolis算法生成一个状态序列，然后
对序列中的每一个状态进行上面的扩充，最终可以生成符合周期条件要求
的f (i)gN
i=1，可以按规则求解能量函数。
关于能量函数的生成，可以使用程序包中的funcH.m函数进行求解。函数
输入格式 （(M 1) (M 1)），输出对应的能量值H。具体的使用方法请参
见程序声明。
3.2 最初的实验结果(Results at the first time)
一切就绪，先来看一组实验结果。运行Ising.m函数（具体参数选取方法见
函数声明）。选取M = 10; N = 100，并采取等可能预选。 取0至1间0.1为间隔
的一组数，分别求出U和C，对 作图。
0 0.2 0.4 0.6 0.8 1
−.5
−
−.5
−
−.5
0
0 0.2 0.4 0.6 0.8 1
−
−.8
−.6
−.4
−.2
0
0.2
0.4
0.6
0.8
1
这样的图形完全让人摸不着头脑。再来看看使用单点变化预选得到的图形。
3
0 0.2 0.4 0.6 0.8 1
−.5
−
−.5
−
−.5
0
0 0.2 0.4 0.6 0.8 1
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
显然，使用单点变化预选的结果要好于等可能预选。
为什么等可能预选的结论会如此糟糕呢？这就要注意到 H的值了。由
于等可能预选会有可能改变很多格点的值，导致 H可能会相当大，这使A =
minf1; exp( H)g 0，从而极大可能会保留原来的 n作为 n+1。这样会生成一
串相同的 列，结论当然糟糕。
为了证明这个论点。取M = 3，此时由于格点数不多，因此，A不会太小，
等可能预选可能可以得到不错的结果。下面就是一组得到的图形组。
0 0.2 0.4 0.6 0.8 1
−
−.5
−
−.5
−
−.5
0
0 0.2 0.4 0.6 0.8 1
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
反观单点变化预选，由于只有一个格点的情况发生了改变，无论M多
大， H的值都不会太大。因此，在后面的讨论中都使用单点变化预选法。
草草的看一下图，可以发现C在0.2到0.4之间有较大的跳跃，对U来讲，也
有部分斜率较大，这有可能成为我们所希望得到的临界点。当然，这个图并不
能完全说明问题。为了找出最能说明结论的图形，我们需要对实验结果进行优
化。
3.3 实验结果的优化(Optimize the Results)
要想让实验结果更加出色，一个最原始但却是十分有效的方法就是多次尝
试。因为Monte-Carlo方法和Metropolis算法都涉及到了随机数，因此每次运行
的结果都不一样，多次尝试可以从中选取比较令人满意的值，达到优化实验的
结果。
当然，不停尝试多少显得有些不上档次，下面讨论一些提高总体图形表现
力手段。它们确实有用，但在此之前，需要先声明的是，由于随机的存在，不
可能出现完美的理想结果，有少数不大支持结论的项是很正常的。
首先，提高N值可以使得Monte-Carlo算法的近似更好，从而起到优化的作
用。举例来看，如果将前面假设中N取为1000的话，有以下的图形，可以说是
个更加令人满意的结果。
4
0 0.2 0.4 0.6 0.8 1
−.5
−
−.5
−
−.5
0
0 0.2 0.4 0.6 0.8 1
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
除了N之外，还有其它的办法优化实验结果。我们的实验都是从全部格点
为一的初始态开始的，在前面的若干步中，由于每次仅能改换一个数字，因此
随机性不强。为避免这样的问题，可以将前N1项去掉，不进入时间序列估计的
范畴。另外，由于相邻两项有很高的关联度，因此，我们可以每隔N2项取一个
状态记录，而不是每项都记录。这样也许会有好处。
下面就以N1 = 800; N2 = 10得到的图像。
0 0.2 0.4 0.6 0.8 1
−.5
−
−.5
−
−.5
0
0.5
0 0.2 0.4 0.6 0.8 1
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
这个图像的效果又有了提升，当然，计算速度更慢，因为原来只需要迭
代到1000项的 就可以了，现在需要800 + 1000 10 = 18000项才能得到最终的
解。看来，一个更好的结果需要付出更大的代价。
我们无法否认的是，无论N1和N2如何取值，各项之间还是有关系的，所
以，要想减小随机偏离，还有一个很基本的方法，就是多个结果求均值。当
然，这又是更费时间的一种策略，但的确是有效果的。同时兼顾时间的复杂
性，我们可以在重复次数和N1; N2 的取值上进行一些取舍和兼顾，生成一组混
合策略。比如，取N = 1000; N1 = 500; N2 = 5; times = 10，可以得到以下图像。
0 0.2 0.4 0.6 0.8 1
−.5
−
−.5
−
−.5
0
0.5
0 0.2 0.4 0.6 0.8 1
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
5
这组图像的中的U的突变更好的显示了出来。
总结一下，要得到一幅好的图像，我们可以通过牺牲时间来达到。采取的
措施包括增大时齐列的长度，除去开头的若干项，隔项选取抽样点，多个结果
求均值等。需要注意的是，无论用什么样的方法，由于具有随机性，每次得到
的结果都不一样，因此，多次尝试并选取较好的图是必要的。
3.4 实验结果的加细(Thinning the Results)
通过上面的实验结果，我们可以大致估计 临界值大概在0.4的周围。但由
于我们仅仅计算了十个不同温度下的内能和比热容的值，无法更进一步的估计
临界值。要想更好的估计，需要对试验结果进行加细，如 取0至1之间以0.02为
间隔的一组数。选取策略N = 1000; N1 = 800; N2 = 10; times = 1进行计算，得到
下面的图形。
0 0.2 0.4 0.6 0.8 1
−.5
−
−.5
−
−.5
0
0.5
0 0.2 0.4 0.6 0.8 1
0
0.2
0.4
0.6
0.8
1
1.2
1.4
加细后的图像变得震荡起来，这是很正常的，因为增加了温度点，这点可
能偏离理想值，从而形成一个尖角。由于温度点的增多，要想通过不断尝试使
每个点都恰好取到理想的地方是相当困难的，因此，我们只有一方面牺牲更多
的时间，一方面不去追求完美，而只是观察大致的走向，得到最终的结果。
取N = 3000; N1 = 2000; N2 = 10; times = 1得到下面的图形。
0 0.2 0.4 0.6 0.8 1
−.5
−
−.5
−
−.5
0
0.5
0 0.2 0.4 0.6 0.8 1
0
0.2
0.4
0.6
0.8
1
1.2
1.4
这组策略运行时间相当长，达到数分钟，但得到的结果已经相当令人满意
了，毕竟，一共有51个随机产生的点啊。
3.5 格点边长的增加(Increase the length of the points)
前面讨论了M = 10的情形，现在看一看扩大M后的结果。方法完全一样，
但是值得注意的是，格点增多后，相应的N; N1; N2都应该增大以达到相同的效
果。当然，这会耗费更长的时间。因此，这里只做出一个较为粗糙的结果，仅
仅作为一种参考。取 间隔0.1，M = 100; N = 1000; N1 = 800; N2 = 5; times = 1。
6
0 0.2 0.4 0.6 0.8 1
−.2
−
−.8
−.6
−.4
−.2
−
−.8
−.6
0 0.2 0.4 0.6 0.8 1
0
1
2
3
4
5
6
7
8
9
10
4 临界点温度(The temperature in the critical point)
4.1 求解策略(Strategy of the solvation)
由于每张图都不大一样，因此，通过图像只能确定一个大概的临界点温
度，比如， = 0:4就是一个大概的点。但由于加细后的震荡变得严重，我们
很难说服自己将一幅图的结论作出最终的结果。因此，我们需要大量的计算结
果。
我们通过下面的策略来求临界点的温度。选取 从0到1以0.02为间隔。选取
参数策略N = 500; N1 = 500; N2 = 3; times = 1。选取得到的51个C值中的最大点
对应的 作为临界点。连续运行50次，取平均值即为临界点温度。
4.2 数值结果(Numerical results)
尽管设定的参数值都不太大，但由于要循环50次，因此速度仍然比较慢。
特别是当M较大时，运行速度相当慢（M = 200是运行了近3小时）。因此，我
们没有足够的时间对所有的M进行实验，仅选取M = 10; 20; 50; 100; 200几个比
较特别的数进行了实验。
下面的组图是使用上面的策略得到的50个 估计值的分布（蓝线）以及平
均的 值（红线）。这可以由solvebeta.m函数导出。由于速度相当的慢，在实
际操作中，根据大致的结果估计， 仅选取了0至0.5中的0.02为间隔的总共26个
值。
M = 10 M = 20
0 5 10 15 20 25 30 35 40 45
0
0.05
0.1
0.15
0.2
0.25
0.3
0.35
0.4
0.45
0.5
Experimental value
Average value
0 5 10 15 20 25 30 35 40 45
0
0.05
0.1
0.15
0.2
0.25
0.3
0.35
0.4
0.45
0.5
Experimental value
Average value
7
M = 50 M = 100
0 5 10 15 20 25 30 35 40 45
0
0.05
0.1
0.15
0.2
0.25
0.3
0.35
0.4
0.45
0.5
Experimental value
Average value
0 5 10 15 20 25 30 35 40 45
0
0.05
0.1
0.15
0.2
0.25
0.3
0.35
0.4
0.45
0.5
Experimental value
Average value
M = 200
0 5 10 15 20 25 30 35 40 45
0
0.05
0.1
0.15
0.2
0.25
0.3
0.35
0.4
0.45
0.5
Experimental value
Average value
得到的 的估计值总结为下表。
M 10 20 50 100 200
M 0.3908 0.3344 0.2004 0.1552 0.1356
从上面的图像和图表结论可以发现，随着M的增大，临界点 的值会发生变
化，并且趋势是逐渐下降。如果认定啊 = (kBT)1 的话，临界点温度是逐渐上
升的。
另外，我们还可以发现，M越大，利用上述策略求得的值的振荡越小，结
果越具有可信度（方差很小）。当M ! 1时，会有一个极限的值 1，这个值是
多少尚不知道，需要进行更多步骤的数值实验才能得到。当然，根据上面的数
据，可以发现 100和 200的值已经相当接近了，可见下降的趋势在减弱，因此可
以大胆估计 1距离 200 不会太远。参考文献[4]中称 1 = 1=8，是有能够让人接
受的。
1反映了全空间范围下的铁磁相变临界点的位置，对于很多物理方面的研
究是相当有帮助的。
5 总结(Summary)
用二维Ising模型模拟铁磁系统的相变是相当成功的。我们得到了很好的图
形结果，并且对于不同边长的系统分别求出了值得信服的临界点值，并进行了
8
对比。如果有更多的时间，还可以进行更多次的数值实验，从而得到关于临界
点的更多有用的性质。
9
A 参考文献(Bibliography)
[1]张平文、李铁军，《数值分析》，2007年1月，北京大学出版社
[2]TEXGuru，《Manual of LATEX2 》，网络中文版
[3]程卫国、冯峰、王雪梅、刘艺，《Matlab5.3精要、编程及高级应用》（网络
版），2000年4月，机械工程出版社
[4]孙前芳，”二维依辛模型相变附近临界指数 的计算机模拟计算”，《科学技
术与工程》，2005.17，p1235-1237
[5]林旭升，”二维依辛模型相变临界点温度的模拟计算”，《大学物理》，2000.5，p13-
15
B 软件支持(Software support)
B.1 数值计算工具(Numerical calculation tools)
本次报告主要使用的数值计算工具为Matlab。通过调用m文件可以进行相
关的计算。函数使用方法可见函数的声明。所有函数在电子版中附上。
B.2 作图工具(Mapping tools)
本次报告的图像是通过Matlab绘制的。图像的最终格式为*.eps，由于此
格式不支持中文导出，故图像中的文字均用了英文批注。部分图形生成后
用Illustrator进行了少量修改。所有图像文件在电子版中附上。
B.3 报告排版工具(Typesetting tools)
本次报告使用了LATEX进行排版，最终电子版以只读型pdf格式上交，同时附
有源文件*.tex供参考和修改。
C 感谢(Thanks)
本次报告能够顺利完成，非常感谢蔡啸天同学提供他的电脑，夜以继日的
运行程序。特别是在求作 100、 200等大型问题时，往往需要超过1小时的连续
运行。感谢他和他的电脑作出的努力，让我得到了一个比较好的结论。
10
Contents
1 介绍(Introduction) 1
2 相变问题与Ising模型
(The problem of phase change and Ising model) 1
2.1 相变问题(The problem of phase change) . . . . . . . . . . . . . . . 1
2.2 Ising模型初探(The first approach to Ising model) . . . . . . . . . . 1
2.3 遍历性(Ergodicity) . . . . . . . . . . . . . . . . . . . . . . . . . . 2
2.4 Metropolis算法(Metropolis algorithm) . . . . . . . . . . . . . . . . 2
2.5 小结(Brief summary) . . . . . . . . . . . . . . . . . . . . . . . . . 3
3 数值实验(Numerical experiment) 3
3.1 生成能量函数(Creating power function) . . . . . . . . . . . . . . 3
3.2 最初的实验结果(Results at the first time) . . . . . . . . . . . . . . 3
3.3 实验结果的优化(Optimize the Results) . . . . . . . . . . . . . . . 4
3.4 实验结果的加细(Thinning the Results) . . . . . . . . . . . . . . . 6
3.5 格点边长的增加(Increase the length of the points) . . . . . . . . . 6
4 临界点温度(The temperature in the critical point) 7
4.1 求解策略(Strategy of the solvation) . . . . . . . . . . . . . . . . . 7
4.2 数值结果(Numerical results) . . . . . . . . . . . . . . . . . . . . . 7
5 总结(Summary) 8
A 参考文献(Bibliography) 10
B 软件支持(Software support) 10
B.1 数值计算工具(Numerical calculation tools) . . . . . . . . . . . . . 10
B.2 作图工具(Mapping tools) . . . . . . . . . . . . . . . . . . . . . . 10
B.3 报告排版工具(Typesetting tools) . . . . . . . . . . . . . . . . . . . 10
C 感谢(Thanks) 10
11
Abhi Pokhrel