计算物理 蒙特卡罗方法基础
合集下载
相关主题
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
实验方案： 在平滑桌面上划一组相距为s的平行线，向此桌面随意地投掷长度l的细
针，那末从针与平行线相交的概率就可以得到π的数值。
SSS L
9
数学统计理论计算：
针的投影长度
Lcos
S
确定的 的概率
，α相交
Lcos
p
cos
s
A
B
0
cos 的平均值
1
cos d
2
M
0
N
假如在N次投针中，有M次和平行线相交。当N充分大时，相交的频数 M/N 就近似为细针与平行线相交的概率。
end do end
17
结果和分析
(1) 总计投点1.0×105次 (2) 该算法收敛，
计算值平均值为3.1392
1.6422/ n 3.1392 1.6422/ n
18
例3 定积分计算
1
I 0 f ( x)dx
0 x1 0 y1
y 1
O
1x
这时我们可以随机地向正方形内投点，最后统计落在曲线下的点数M， 当总的掷点数N充分大时，M/N就近似等于积分值I。
这意味着试验所得的值的不确定性的范围如下：
对100次投针为，
0.2374
对10,000次投针为， 0.0237
对1,000,000次投针为， 0.0024
可见，增加模拟的次数可以减小误差，但不可消除误差。
12
前人进行了实验，其结果列于下表 ：
实验者
年份
投计次数
π的实验值
沃尔弗(Wolf)
1850
5000
例如：有一个随机变量η，它满足分布密度函数f(x)。如果我们将n个满足分布 密度函数f(x)的独立随机数相加：
Rn 1 2 n
则Rn满足高斯分布。高斯分布可以由给定的期望值μ和方差σ完全确定下来。
N( , 2 ) 1 exp[ ( x )2 / 2 2 ] 2
28
lim
n
1 n
(2) x , y 为独立变量
协方差=0
V{x y} V{x}V{ y}
26
五 大数法则和中心极限定理
概率论中的大数法则和中心极限定理是蒙特卡洛方法的基础。
1 大数法则反映了大量随机数之和的性质。
如果函数在[a,b]区间，以均匀的概率分布密度随机地取n个数ui，对每个计 算出函数值h[ui]。大数法则告诉我们这些函数值之和除以n所得的值将收敛于函 数h在区间[a,b]的期望值，即
的标准误差 4 p(1 p) / n 4 0.4155/ n 1.6422/ n
对100次投针为，
0.1642
对10,000次投针为， 0.0164
对1,000,000次投针为， 0.0016
15
投点法实验程序流程图
n n max
Yes
n n1
产生随机数
1 ,2
x L1 , y L2
gN
1 N
N
g(ri )
i 1
代表了该运动员的成绩。换言之， <g>为积分的估计值，或近似值。
假设射击100次，
环数
7
8
9 10
击中次数 10 10 30 50
概率
0.1 0.1 0.3 0.5
平均成绩
gN
1 N
N
920
i1 g(ri ) 100
9.2
7
设r表示射击运动员的弹着点到靶心的距离，g(r)表示击中r处相应的得分 数（环数），f(r)为该运动员的弹着点的分布密度函数，它反映运动员的射 击水平。该运动员的射击成绩为
这种模拟是以所谓“马尔科夫”（Markov）链的形式产生系统的分 布序列。该方法可以使我们能够研究经典和量子多粒子系统的问题。
5
一 基本思想
直接蒙特卡洛模拟法： 对求解问题本身就具有概率和统计性的情况。
如：中子在介质中的传播，核衰变过程等， 思想是按照实际问题所遵循的概率统计规律，用计算机进行直接的抽样
(2)分子动力学方法：确定性模拟方法。它是通过数值求解一个个
的粒子运动方程来模拟整个系统的行为。在统计物理中称为分子动力学 （Molecular Dynamics）方法。
(3)离散型模拟方法 -- 元胞自动机等
1
2-1 蒙特卡罗方法的基础知识
What is a Monte Carlo method ?
Metropolis coined the name “Monte Carlo”, from its gambling casino.
4
从蒙特卡洛模拟的应用来看，该类型的应用可以分为三种形式： （1）直接蒙特卡洛模拟。
它采用随机数序列来模拟复杂随机过程的效应。 （2）蒙特卡洛积分。
这是利用随机数序列计算积分的方法。积分维数越高，该方法的积分 效率就越高。 （3）Metropolis蒙特卡洛模拟
24
一个函数或变量的方差:
V { f } E{( f E{ f })2 }
标准误差或均方根误差：方差的平方根。
V( f )
由于标准误差与其真值有相同的量纲，因而它比方差更具有物理意义。
假如x和y是随机变量，c是一个常数，则：
E{cx y} cE{x} E{ y} V {cx y} c2V { x} V { y}
(1)数学期望是线性算符 (2)方差是非线性算符
ห้องสมุดไป่ตู้
2cE{(x E{x})(y E{ y})}
x , y间的协方差
25
V {cx y} c2V { x} V { y} 2cE{(x E{x})(y E{ y})}
协方差>0，正关联 协方差<0，负关联
注意：
(1) 协方差=0
x , y 为独立变量
试验，然后计算其统计参数。 该方法也就是通常所说的“计算机实验”。
间接蒙特卡洛方法： 蒙特卡洛方法也可以人为地构造出一个合适的概率模型，依照该模型进
行大量的统计实验，使它的某些统计参量正好是待求问题的解。
6
例1 射击问题（打靶游戏）-- 直接蒙特卡洛方法
现假设该运动员进行了N次射击，每次射击的弹着点依次为r1，r2，…，rN， 则N次得分g(r1)，g(r2)，…，g(rN)的算术平均值
the Comte de Buffon needle experiment, AD 1777 SSS L
2
蒙特卡洛方法的起源
Stanislaw Ulam (1909-1984)
Nicholas Metropolis (1915-1999)
3
The Name of the Game
Monte-Carlo, Monaco
n i 1
h(ui )
lim
n
I
n
1 ba
b
h(u)du I
a
N( , 2 ) 1 exp[ ( x )2 / 2 2 ] 2
当n充分大时，对任意的λ，由列维定理知：
lim
n
P{
(f n
g 0 g(r) f (r)dr
用概率语言来说， g(r) 是随机变量， <g>的数学期望，即
g Eg(r) 7*0.1 8*0.1 9*0.3 10*0.5
在该例中，用N次试验所得成绩的算术平均值作为数学期望<g>的估计值 （积分近似值）。
8
例2 圆周率的数值计算 -- 间接蒙特卡洛方法 (1) 巴夫昂(Buffon)投针实验
2N
M
10
经过n次投针后得到π值的精度
针与平行线相交的概率
p 2
针与平行线相交的次数应满足二项式分布
相交和不相交
其期望值为
np
2 的方差
2 p(1 p) / n
2 的标准误差
的标准误差
p(1 p) / n
2
2
p(1 p) / n 2.37 / n
11
的标准误差
2
2
p(1 p) / n 2.37 / n
投点法实验源程序
do i=1,nmax call randomnum() x=lenr*(rand-0.5d0)*2.0d0 call randomnum() y=lenr*(rand-0.5d0)*2.0d0 dxy2=x*x+y*y if ( dxy2 .le. lenr2 ) then m=m+1 end if ncount=ncount+1 if (mod(ncount,100) .eq. 0 ) then write(10,"(I10,F15.6)")ncount, 4.0d0*dble(m)/dble(ncount) end if
2L 任意投点落在圆内的概率
p Scircle S square
L2
(2L)2
4
14
投针实验的误差分析
2 的标准误差
的标准误差
p(1 p) / n
2
2
p(1 p) / n 2.37 / n
投点实验的误差分析
的标准误差 p(1 p) / n
4
p Scircle
Ssquare 4
19
直接蒙特卡罗方法的思想
该方法是按照实际问题所遵循的概率统计规律，用计算机进行直接的抽 样试验，然后计算其统计参数。
间接蒙特卡罗方法的思想s
当问题可以抽象为某个确定的数学问题时， (1)首先建立一个恰当的概率模型，即确定某个随机事件A或
随机变量X，使得待求的解等于随机事件出现的概率或随 机变量的数学期望值。 (2)然后进行模拟实验，即重复多次地模拟随机事件A或随机 变量X。 (3)最后对随机实验结果进行统计平均，求出A出现的频数或 X的平均值作为问题的近似解。
3.1596
斯密思(Smith)
1855
3204
3.1553
福克斯(Fox)
1894
1120
3.1419
拉查里尼(Lazzarini) 1901
3408
3.1415929
13
(2) 投点法实验
实验方案： 在平滑桌面上划正方形，同时划一内切圆，向此正方形随意地投点，那
末投点落在圆内的概率就可以得到的π数值。
1 b
1n
I ba
a
h(u)du
lim
n
n
i 1
h(ui )
lim
n
I
n
大数法则保证了在抽取足够多的随机样本后，计算得到的积分的蒙特卡洛估 计值将收敛于该积分的正确结果。
27
2 中心极限定理
中心极限定理告诉我们：在有足够大，但又有限的抽样数n的情况下，蒙特卡 洛估计值是如何分布的。
该定理指出：无论随机变量的分布如何，它的若干个独立随机变量抽样值之和 总是满足正则分布(即高斯分布)。
计算机模拟方法
(1)蒙特卡洛方法：随机性模拟方法或统计试验方法，又称蒙特卡洛
(Monte Carlo)方法。它是通过不断产生随机数序列来模拟过程。自然界中有 的过程本身就是随机的过程，物理现象中如粒子的衰变过程、粒子在介质中的 输运过程...等。当然蒙特卡洛方法也可以借助概率模型来解决不直接具有随机 性的确定性问题。
若h(u,v)=p(u) ·q(v)，则两个随机变量u’和v’彼此独立。
对如下三个变量
xr ys z mod(( r s),1)
(x , y)彼此独立； (x , z)彼此独立； (y , z)彼此独立；
(x , y , z)相互关联。
23
四 期望值、方差和协方差
一个函数f(u’) 的数学期望值定义为该函数的平均值
u’在u到u+du之间值的概率。
G(u) 称为g(u)的分布函数。
x
G(u) g( x)dx
G(u)在区间取值的单调递增函数
g(u) dG(u) / du
b
a g( x)dx 1
22
三 随机变量的独立性
假如我们考虑两个随机变量u’和v’的分布，则必须引进这两个变量的联合 分布密度函数h(u,v)，此时带来的数学问题就更为复杂。
E{ f } f (u)dG(u) f (u)g(u)du
G(u) 称为u的分布函数。
通常u是在[a,b]区间均匀分布的随机变量，有
dG(u) du /(b a)
f 的数学期望值：
E{
f
}
b
1 a
f
(u)du
类似地，自由变量u的期望值为u的平均值，有
E{u} udG(u) ug(u)du
r 2 x2 y2 L2
Yes
M M1
计 算
4M dble( )
N
16
program main use freconstant use randomname implicit none integer nmax,m integer i,ncount real*8 lenr,lens,lenr2 real*8 x,y,dxy2 open(10,file='Pi.dat') call randomval() lenr=1.0d0 lens=2.0d0*lenr lenr2=lenr*lenr m=0 ncount=0 write(*,*)"Input nmax:" read(*,*)nmax
20
作业
“Buffon投针法”计算圆周率。
21
二 随机变量和随机变量的分布
随机变量：是一个不止是一个值的变量（通常是连续的）， 并且人们可能无法事先预言某一个特定的值。
但是：其分布是可以了解的，假设我们研究某一连续性的变量，由随机变量 的分布我可以得到它取某值的概率：
P(u u' u du) g(u)du g(u) 称为u的概率分布密度函数，它表示随机变量
针，那末从针与平行线相交的概率就可以得到π的数值。
SSS L
9
数学统计理论计算：
针的投影长度
Lcos
S
确定的 的概率
，α相交
Lcos
p
cos
s
A
B
0
cos 的平均值
1
cos d
2
M
0
N
假如在N次投针中，有M次和平行线相交。当N充分大时，相交的频数 M/N 就近似为细针与平行线相交的概率。
end do end
17
结果和分析
(1) 总计投点1.0×105次 (2) 该算法收敛，
计算值平均值为3.1392
1.6422/ n 3.1392 1.6422/ n
18
例3 定积分计算
1
I 0 f ( x)dx
0 x1 0 y1
y 1
O
1x
这时我们可以随机地向正方形内投点，最后统计落在曲线下的点数M， 当总的掷点数N充分大时，M/N就近似等于积分值I。
这意味着试验所得的值的不确定性的范围如下：
对100次投针为，
0.2374
对10,000次投针为， 0.0237
对1,000,000次投针为， 0.0024
可见，增加模拟的次数可以减小误差，但不可消除误差。
12
前人进行了实验，其结果列于下表 ：
实验者
年份
投计次数
π的实验值
沃尔弗(Wolf)
1850
5000
例如：有一个随机变量η，它满足分布密度函数f(x)。如果我们将n个满足分布 密度函数f(x)的独立随机数相加：
Rn 1 2 n
则Rn满足高斯分布。高斯分布可以由给定的期望值μ和方差σ完全确定下来。
N( , 2 ) 1 exp[ ( x )2 / 2 2 ] 2
28
lim
n
1 n
(2) x , y 为独立变量
协方差=0
V{x y} V{x}V{ y}
26
五 大数法则和中心极限定理
概率论中的大数法则和中心极限定理是蒙特卡洛方法的基础。
1 大数法则反映了大量随机数之和的性质。
如果函数在[a,b]区间，以均匀的概率分布密度随机地取n个数ui，对每个计 算出函数值h[ui]。大数法则告诉我们这些函数值之和除以n所得的值将收敛于函 数h在区间[a,b]的期望值，即
的标准误差 4 p(1 p) / n 4 0.4155/ n 1.6422/ n
对100次投针为，
0.1642
对10,000次投针为， 0.0164
对1,000,000次投针为， 0.0016
15
投点法实验程序流程图
n n max
Yes
n n1
产生随机数
1 ,2
x L1 , y L2
gN
1 N
N
g(ri )
i 1
代表了该运动员的成绩。换言之， <g>为积分的估计值，或近似值。
假设射击100次，
环数
7
8
9 10
击中次数 10 10 30 50
概率
0.1 0.1 0.3 0.5
平均成绩
gN
1 N
N
920
i1 g(ri ) 100
9.2
7
设r表示射击运动员的弹着点到靶心的距离，g(r)表示击中r处相应的得分 数（环数），f(r)为该运动员的弹着点的分布密度函数，它反映运动员的射 击水平。该运动员的射击成绩为
这种模拟是以所谓“马尔科夫”（Markov）链的形式产生系统的分 布序列。该方法可以使我们能够研究经典和量子多粒子系统的问题。
5
一 基本思想
直接蒙特卡洛模拟法： 对求解问题本身就具有概率和统计性的情况。
如：中子在介质中的传播，核衰变过程等， 思想是按照实际问题所遵循的概率统计规律，用计算机进行直接的抽样
(2)分子动力学方法：确定性模拟方法。它是通过数值求解一个个
的粒子运动方程来模拟整个系统的行为。在统计物理中称为分子动力学 （Molecular Dynamics）方法。
(3)离散型模拟方法 -- 元胞自动机等
1
2-1 蒙特卡罗方法的基础知识
What is a Monte Carlo method ?
Metropolis coined the name “Monte Carlo”, from its gambling casino.
4
从蒙特卡洛模拟的应用来看，该类型的应用可以分为三种形式： （1）直接蒙特卡洛模拟。
它采用随机数序列来模拟复杂随机过程的效应。 （2）蒙特卡洛积分。
这是利用随机数序列计算积分的方法。积分维数越高，该方法的积分 效率就越高。 （3）Metropolis蒙特卡洛模拟
24
一个函数或变量的方差:
V { f } E{( f E{ f })2 }
标准误差或均方根误差：方差的平方根。
V( f )
由于标准误差与其真值有相同的量纲，因而它比方差更具有物理意义。
假如x和y是随机变量，c是一个常数，则：
E{cx y} cE{x} E{ y} V {cx y} c2V { x} V { y}
(1)数学期望是线性算符 (2)方差是非线性算符
ห้องสมุดไป่ตู้
2cE{(x E{x})(y E{ y})}
x , y间的协方差
25
V {cx y} c2V { x} V { y} 2cE{(x E{x})(y E{ y})}
协方差>0，正关联 协方差<0，负关联
注意：
(1) 协方差=0
x , y 为独立变量
试验，然后计算其统计参数。 该方法也就是通常所说的“计算机实验”。
间接蒙特卡洛方法： 蒙特卡洛方法也可以人为地构造出一个合适的概率模型，依照该模型进
行大量的统计实验，使它的某些统计参量正好是待求问题的解。
6
例1 射击问题（打靶游戏）-- 直接蒙特卡洛方法
现假设该运动员进行了N次射击，每次射击的弹着点依次为r1，r2，…，rN， 则N次得分g(r1)，g(r2)，…，g(rN)的算术平均值
the Comte de Buffon needle experiment, AD 1777 SSS L
2
蒙特卡洛方法的起源
Stanislaw Ulam (1909-1984)
Nicholas Metropolis (1915-1999)
3
The Name of the Game
Monte-Carlo, Monaco
n i 1
h(ui )
lim
n
I
n
1 ba
b
h(u)du I
a
N( , 2 ) 1 exp[ ( x )2 / 2 2 ] 2
当n充分大时，对任意的λ，由列维定理知：
lim
n
P{
(f n
g 0 g(r) f (r)dr
用概率语言来说， g(r) 是随机变量， <g>的数学期望，即
g Eg(r) 7*0.1 8*0.1 9*0.3 10*0.5
在该例中，用N次试验所得成绩的算术平均值作为数学期望<g>的估计值 （积分近似值）。
8
例2 圆周率的数值计算 -- 间接蒙特卡洛方法 (1) 巴夫昂(Buffon)投针实验
2N
M
10
经过n次投针后得到π值的精度
针与平行线相交的概率
p 2
针与平行线相交的次数应满足二项式分布
相交和不相交
其期望值为
np
2 的方差
2 p(1 p) / n
2 的标准误差
的标准误差
p(1 p) / n
2
2
p(1 p) / n 2.37 / n
11
的标准误差
2
2
p(1 p) / n 2.37 / n
投点法实验源程序
do i=1,nmax call randomnum() x=lenr*(rand-0.5d0)*2.0d0 call randomnum() y=lenr*(rand-0.5d0)*2.0d0 dxy2=x*x+y*y if ( dxy2 .le. lenr2 ) then m=m+1 end if ncount=ncount+1 if (mod(ncount,100) .eq. 0 ) then write(10,"(I10,F15.6)")ncount, 4.0d0*dble(m)/dble(ncount) end if
2L 任意投点落在圆内的概率
p Scircle S square
L2
(2L)2
4
14
投针实验的误差分析
2 的标准误差
的标准误差
p(1 p) / n
2
2
p(1 p) / n 2.37 / n
投点实验的误差分析
的标准误差 p(1 p) / n
4
p Scircle
Ssquare 4
19
直接蒙特卡罗方法的思想
该方法是按照实际问题所遵循的概率统计规律，用计算机进行直接的抽 样试验，然后计算其统计参数。
间接蒙特卡罗方法的思想s
当问题可以抽象为某个确定的数学问题时， (1)首先建立一个恰当的概率模型，即确定某个随机事件A或
随机变量X，使得待求的解等于随机事件出现的概率或随 机变量的数学期望值。 (2)然后进行模拟实验，即重复多次地模拟随机事件A或随机 变量X。 (3)最后对随机实验结果进行统计平均，求出A出现的频数或 X的平均值作为问题的近似解。
3.1596
斯密思(Smith)
1855
3204
3.1553
福克斯(Fox)
1894
1120
3.1419
拉查里尼(Lazzarini) 1901
3408
3.1415929
13
(2) 投点法实验
实验方案： 在平滑桌面上划正方形，同时划一内切圆，向此正方形随意地投点，那
末投点落在圆内的概率就可以得到的π数值。
1 b
1n
I ba
a
h(u)du
lim
n
n
i 1
h(ui )
lim
n
I
n
大数法则保证了在抽取足够多的随机样本后，计算得到的积分的蒙特卡洛估 计值将收敛于该积分的正确结果。
27
2 中心极限定理
中心极限定理告诉我们：在有足够大，但又有限的抽样数n的情况下，蒙特卡 洛估计值是如何分布的。
该定理指出：无论随机变量的分布如何，它的若干个独立随机变量抽样值之和 总是满足正则分布(即高斯分布)。
计算机模拟方法
(1)蒙特卡洛方法：随机性模拟方法或统计试验方法，又称蒙特卡洛
(Monte Carlo)方法。它是通过不断产生随机数序列来模拟过程。自然界中有 的过程本身就是随机的过程，物理现象中如粒子的衰变过程、粒子在介质中的 输运过程...等。当然蒙特卡洛方法也可以借助概率模型来解决不直接具有随机 性的确定性问题。
若h(u,v)=p(u) ·q(v)，则两个随机变量u’和v’彼此独立。
对如下三个变量
xr ys z mod(( r s),1)
(x , y)彼此独立； (x , z)彼此独立； (y , z)彼此独立；
(x , y , z)相互关联。
23
四 期望值、方差和协方差
一个函数f(u’) 的数学期望值定义为该函数的平均值
u’在u到u+du之间值的概率。
G(u) 称为g(u)的分布函数。
x
G(u) g( x)dx
G(u)在区间取值的单调递增函数
g(u) dG(u) / du
b
a g( x)dx 1
22
三 随机变量的独立性
假如我们考虑两个随机变量u’和v’的分布，则必须引进这两个变量的联合 分布密度函数h(u,v)，此时带来的数学问题就更为复杂。
E{ f } f (u)dG(u) f (u)g(u)du
G(u) 称为u的分布函数。
通常u是在[a,b]区间均匀分布的随机变量，有
dG(u) du /(b a)
f 的数学期望值：
E{
f
}
b
1 a
f
(u)du
类似地，自由变量u的期望值为u的平均值，有
E{u} udG(u) ug(u)du
r 2 x2 y2 L2
Yes
M M1
计 算
4M dble( )
N
16
program main use freconstant use randomname implicit none integer nmax,m integer i,ncount real*8 lenr,lens,lenr2 real*8 x,y,dxy2 open(10,file='Pi.dat') call randomval() lenr=1.0d0 lens=2.0d0*lenr lenr2=lenr*lenr m=0 ncount=0 write(*,*)"Input nmax:" read(*,*)nmax
20
作业
“Buffon投针法”计算圆周率。
21
二 随机变量和随机变量的分布
随机变量：是一个不止是一个值的变量（通常是连续的）， 并且人们可能无法事先预言某一个特定的值。
但是：其分布是可以了解的，假设我们研究某一连续性的变量，由随机变量 的分布我可以得到它取某值的概率：
P(u u' u du) g(u)du g(u) 称为u的概率分布密度函数，它表示随机变量