ising01 能量函数，满足Gibbs分布 相邻项的乘积之和

数值代数上机作业第七次

二维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( )i􀀀hH( )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(M􀀀1)

:::

: : :

:::

(M􀀀1)1 (M􀀀1)(M􀀀1)

1CCCCCCCCCCA

=)

0BBBBBBBBBBBBBBB@

11 1(M􀀀1) 11

:::

: : :

:::

:::

(M􀀀1)1 (M􀀀1)(M􀀀1) (M􀀀1)1

11 1(M􀀀1) 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