计算物理第一章讲义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,02d x ϕπ≤≤≤≤. 而针与一直线相交的充要条件是：sin 2l 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)：cleard=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<=nhold onxi=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 onh_point2=plot(xx,yy,'.r','LineWidth',1); pause(0.01)%暂停0.5秒set(h_point1,'color','k')%将当前直线的颜色由原来的颜色变为黑色k=k+1;frame(k)=getframe(gcf);endwritegif(' 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==1IMWRITE(im,map2,name,'GIF',…'overwrite','DelayTime',dt,'LoopCount',inf);elseIMWRITE(im,map2,name,…'WriteMode','append',…'DelayTime',dt); %,'LoopCount',inf);end end设2,1d l ==,(总投针数n=10000000). Matlab 实际计算程序(CH1.1.2)： cleard=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;y=0.5*l*sin(xi); %产生曲线上对应xi的函数值y for ii =1: nif yi(ii)<=y(ii)m =m+1; %统计处在曲边梯形内的随机点数目endendPi=n/m共计算10次, 其结果见表1.1.1.表1.1.1计算次数 1 2 3 4 5计算值 3.1443 3.1419 3.1427 3.1441 3.1426计算次数 6 7 8 9 10计算值 3.1428 3.1399 3.1430 3.1425 3.14141.1.2计算圆周率的其它方法著名的投针问题利用单位圆与边长为1的正方形面积之比来计算 的近似值. 具体思想如下:如图 1.1.4所示, 单位圆的14为一个扇形G, 它是边长为1的正方形的一部分. 考虑扇形面积在正方形π, 然后乘以4面积中所占的比例k, 得出其结果为4就可以得到π的值.这里如何计算比例k, 运用蒙特卡罗方法的随机投点思想. 在正方形中随机投入很多点, 使所投点几何概率在正方形中每一个位置的机会均等, 然后考察有多少点落在扇形内. 其中落在扇形内的点的个数m与投点总数n之比就是k的近似值.模拟实验程序(CH1.1.3)：cleard=1;n=1000;%绘制几何图形x=linspace(0,1,1000);y=sqrt(1-x.^2);plot(x,y,'-r','LineWidth',2);hold onx1=linspace(0,1,500);y1=1;x2=1;y2=linspace(0,1,500);plot(x1,y1,'-k',x2,y2,'-k','LineWidth',2) axis([0 1 0 1])set(gca,'XTick',0: 1)hold on%模拟投掷实验k=0;while k<=nhold onxi=unifrnd(0,1,1,1);yi=unifrnd(0,1,1,1);plot(xi,yi,'.r')h_point=plot(xi,yi,'.r','LineWidth',1); pause(0.0001) %暂停0.0001秒set(h_point,'color','b')%将当前点的颜色由原来的颜色变为蓝色k=k+1;frame(k)=getframe(gcf);endwritegif(' Throwing experiment.gif',frame,0.1); %计算pi值xi=unifrnd(0,1,n,1);yi=unifrnd(0,1,n,1);%产生n个符合题意随机点的坐标(xi,yi)m=0;for ii =1: nif yi(ii)^2+xi(ii)^2<=1m=m+1;%统计处在曲边梯形内的随机点数目endendmm=mnn=nPi=m/n*4所调用的函数程序如浦丰投针实验.Matlab实际计算程序(CH1.1.4)：cleard=1;n=10000000;xi=unifrnd(0,1,n,1);yi=unifrnd(0,1,n,1);%产生n个符合题意随机点的坐标(xi,yi)k=0;for ii =1: nif yi(ii)^2+xi(ii)^2<=1k =k+1; %统计处在曲边梯形内的随机点数目endendPi=k/n*4实际计算投掷次数取n=10000000, 计算结果见表1.1.2.表1.1.2计算次数 1 2 3 4 5计算值 3.1421 3.1415 3.1417 3.1419 3.1423计算次数 6 7 8 9 10计算值 3.1416 3.1410 3.1410 3.1412 3.1416 从上述结果来看, 计算的结果涨落比较大, 如果要达到比较高的精度, 需要大幅增加投掷的次数.1.2利用蒙特卡罗法求解数值积分和函数极值1.2.1定积分的蒙特卡罗算法假定函数()f x 在[],a b 内有界连续且()0f x ≥, 计算定积分()baI f x dx =⎰.为计算出定积分值, 可构造一概率模型: 取一个边长分别为b a -和h 的矩形D (如图 1.2.1), 使曲边梯形在矩形域之内, 并在矩形内随机投点, 假设随机点均匀地落在整个矩形之内, 则落在图中灰色区域内的随机点数k 与投掷点总数N 之比k N , 就近似地等于灰色区面积与矩形面积之比, 从而得出定积分()kI b a h N=-.例1.2.1: 计算211x I edx --=⎰.对()2x f x e-=, 我们很难求出其原函数, 所以用牛顿-莱布尼茨公式无法求解(见图1.2.2).这个问题可运用蒙特卡罗方法依如下步骤求出其近似值.步骤1: 产生矩形域D内的N个均匀随机点Pi(xi,yi);步骤2: 统计满足条件yi≤f(x)的落在曲边梯形内随机点Pi的数目k;步骤3: (取h =1,a =-1, b =1,N =100000), 则定积分的近似值I≈(b-a)hk/N =k/5000.Matlab计算程序(CH1.2.1)：clearxi =unifrnd(-1,1,100000,1);yi =rand(100000,1);%产生100000个符合题意随机点的坐标(xi,yi)k =0;y = exp(-xi.^2); %产生曲线上对应xi的函数值y for i =1: 100000if yi(i) <=y(i)k =k+1; %统计处在曲边梯形内的随机点数目endendI=k/50000将以上程序利用matlab软件重复执行6次, 所得计算结果见表1.2.1.表1.2.1第i次 1 2 3 4 5 6积分近似值1.4934 1.4928 1.4963 1.4936 1.4946 1.4941 注意: 由于随机误差的影响, 每次计算机模拟的结果可能不一样, 但是各次计算结果总是在准确值附近作微小摆动.Jinzhong University1.2.2二重积分的蒙特卡罗算法对二重积分(),I f x y dxdy =⎰⎰, 设(),f x y 为区域A 上的有界函数且(),0f x y ≥，据其几何意义,它是以(),f x y 为曲面顶, A 为底的曲顶柱体C 的体积. 据此, 用均匀随机数计算二重积分的蒙特卡罗方法基本思路为:假设曲顶柱体C 包含在己知体积为D V 的几何体D 的内部, 在D 内产生N 个均匀随机点, 统计出在C 内部的随机点数目N , 则()D N C I V N N ≈.Jinzhong University例1.2.2: 计算(222211A x y x y dxdy ⎡+---+⎢⎣⎰⎰, 其中(){}22,1A A x y x y =+≤. 利用Matlab 先做出函数的图像(见图1.2.3).图 1.2.3分析:该二重积分可看作以222(11)z x x y =+---+为顶的曲顶柱体的体积, 此曲顶柱体在一个边长为2的立方体之内(见图1.2.3), 可计算出其精确值为π.现利用蒙特卡罗算法计算其近似值，见以下程序.Matlab计算程序(CH1.2.2)：clearn =400000;x =unifrnd(-1,1,n,1);y =unifrnd(-1,1,n,1);zi =unifrnd(0,2,n,1);z =1+sqrt(1-x.^2-y.^2) -sqrt(x.^2+y.^2); Nc =0;for ii =1: nif x(ii)^2+y(ii)^2<=1&zi(ii) <=z(ii)Nc =Nc+1;endendI =8*Nc/n下面是重复10次运算所得数据(见表1.2.2)表1.2.2第k次 1 2 3 4 5近似值 3.1429 3.1426 3.1419 3.1453 3.1410第k次 6 7 8 9 10近似值 3.1413 3.1428 3.1393 3.1392 3.1438从上面十个数据可以看出, 每次计算结果在精确值π的附近摆动. 为提高精度, 可选取其它随机数促进蒙特卡罗的收敛性, 例如选取有利随机数来代替均匀随机数可在较少随机点的情况下使计算精度大大提高.1.2.3蒙特卡罗求一元函数的极值高等数学中求一元函数()f x 在闭区间[],a b 的极值的方法为: 找出函数()f x 的所有驻点和不可导点并计算出其相应函数值, 并将其与区间端点函数值比较得出函数的最大值最小值. 但实际问题中方程()0f x '=很难求出其根, 所以驻点往往很难求出. 若利用蒙特卡罗方法: 在[],a b 上产生n 个随机数, 计算出这些数的函数值并作比较, 则可得出极值近似值.例1.2.3: ()()31sin3,22f x x x x ππ=--≤≤, 求函数的最大值.绘制函数图形(图1.2.4).解法一: 使用蒙特卡罗方法Matlab 计算程序(CH1.2.3)：Jinzhong University clearn =1e4;x =unifrnd(-2*pi,2*pi,n,1);fx = (1-x.^3).*sin(3*x);fmax =max(fx)将上述程序重复运行6次(见表1.2.3), 则可得出函数的最大值为194.9062表1.2.3第i次 1 2 3 4 5 6最大值194.906 194.9062 194.9055 194.9056 194.9061 194.9062解法二: 采用等距法搜索函数法求极值Jinzhong University Matlab计算程序(CH1.2.4)：clearn =1e4;x = linspace(-2*pi,2*pi,n);fx = (1-x.^3).*sin(3*x);fmax =max(fx)采用等距法搜索函数极值得函数最大值为194.9061. 根据图 1.2.4, 可知两种方法的结果都很接近函数最大值. 如果取n=1e5所得结果亦为194.9062.。