计算物理第一章讲义PPT

Jinzhong University 计算物理

晋中学院宫建平

Jinzhong University 第1章蒙特卡罗方法的应用

1.1 蒙特卡罗方法简介

1.2 利用蒙特卡罗法求解数值积分和函数极值

*1.3 基于蒙特卡罗的电子双缝衍射的计算机模拟

1.1蒙特卡罗方法简介

蒙特卡罗方法也称随机模拟方法, 有时也称作随机抽样技术或统计实验方法. 它的基本思想是：首先建立一个概率模型或随机过程, 使它的参数等于问题的解; 然后利用计算机模拟该随机现象, 通过对大量模拟仿真试验的结果来分析计算所求参数, 得出实际问题的近似解.

蒙特卡罗方法的特点可归纳成三个方面:

(1)蒙特卡罗方法及其程序结构简单.

(2)蒙特卡罗方法的收敛性及收敛速度与问题的维数无关.

(3) 蒙特卡罗方法的适用性强, 可用在求很多解析方法或常规数值方法难解问题的低精度解.

1.1.1蒲丰投针问题

著名的投针问题是几何概率一个早期的例子, 它是由法国科学家蒲丰(Buffon)在1777年提出的, 因而被称之为蒲丰投针问题. 蒲丰投针问题的解决不仅较典型的反映了几何概率的特征及处理方法, 而且还可以由此了解蒙特卡洛(Monte一Carlo)方法.

蒲丰投针问题: 平面上画有等距离的平行线, 每两条平行线之间的距离为d , 向平面任意投掷一枚长为()l l d <的针, 试求针与平行线相交的概率.

解：设x 表示针落下后针的中点M 到最近的一条平行线的距离, ϕ表示针与平行线所成的角(见图1.1.1), 则

0,02

d x ϕπ≤≤≤≤. 而针与一直线相交的充要条件是：

sin 2

l x ϕ≤.

我们把ϕ及x 表为平面上一点的直角坐标, 则所有基本事件可以用边长为π及2d 的矩形内的点表示出来, 而“针与直线相交”这一事件所包含的基本事件可以用图1.1.2中阴影部分内的点表示出来, 而所求概率

2sin 22l l d P d d πϕϕππ==⎰.

(1) 取白纸一张, 在上面画许多间距为d的等距平行线；

(2)取一根长度为()

l l d

<的均匀直针, 随机地向画有平行线的纸投去, 共投n次, 观察针和直线相交的次数m；

(3)直线与针相交概率P的近似值可用m n得到, 其中n是总投针次数而m是与任一直线相交的总次数, 进而由可得π的近似值为2nl md.

浦丰投针实验模拟程序(CH1.1.1)：

clear

d=2; %平行直线间隔

l=1; %针长

k=0;

n=100; %投掷次数

x1=0: 0.01: 20;

y1=2;y2=4;y3=6;y4=8;y5=10;y6=12;y7=14;y8=16;y 9=18;

plot(x1,y1,'-b',x1,y2,'-b',x1,y3,'-b',x1,y4,'-b',x1,y5,'-

b',x1,y6,'-b',x1,y7,'-b',x1,y8,'-b',x1,y9,'-b','LineWidth', 2)

while k<=n

hold on

xi=unifrnd(2,18,1,1); %生成2和18之间的1 1随机数.

yi=unifrnd(2,18,1,1);

theta=unifrnd(0,2*pi,1,1);

xj=xi+cos(theta);

yj=yi+sin(theta);

xx=xi+0.5*cos(theta);

yy=yi+0.5*sin(theta);

x=[xi,xj];

y=[yi,yj];

zz=0.5*l*sin(theta);

axis([0 20 0 20])

h_point1=plot(x,y,'-g','LineWidth',1);

hold on

h_point2=plot(xx,yy,'.r','LineWidth',1); pause(0.01)%暂停0.5秒

set(h_point1,'color','k')

%将当前直线的颜色由原来的颜色变为黑色k=k+1;

frame(k)=getframe(gcf);

end

writegif(' Needle cast experiment.gif',frame,0.1); 所调用的函数程序如下:

function res=writegif(name,frames,dt)

nframe = length(frames);

for ii=1:nframe

[image,map] = frame2im(frames(ii));

[im,map2]=rgb2ind(image,128);

if ii==1

IMWRITE(im,map2,name,'GIF',…

'overwrite','DelayTime',dt,'LoopCount',inf);

else

IMWRITE(im,map2,name,…

'WriteMode','append',…

'DelayTime',dt); %,'LoopCount',inf);

end end

设2,1d l ==,(总投针数n=10000000). Matlab 实际计算程序(CH1.1.2)： clear

d=2;

l=1;

theta=pi;

n=10000000;

xi=unifrnd(0,theta,n,1);

yi=unifrnd(0,d/2,n,1);

%产生n个符合题意随机点的坐标(xi,yi) m=0;