计算物理第一章讲义

1晋中学院JinzhongUniversity晋中学院JinzhongUniversity晋中学院JinzhongUniversity晋中学院JinzhongUniversity计算物理晋中学院宫建平2晋中学院JinzhongUniversity晋中学院JinzhongUniversity第1章蒙特卡罗方法的应用1.1蒙特卡罗方法简介1.2利用蒙特卡罗法求解数值积分和函数极值*1.3基于蒙特卡罗的电子双缝衍射的计算机模拟3晋中学院JinzhongUniversity晋中学院JinzhongUniversity1.1蒙特卡罗方法简介蒙特卡罗方法也称随机模拟方法,有时也称作随机抽样技术或统计实验方法.它的基本思想是：首先建立一个概率模型或随机过程,使它的参数等于问题的解;然后利用计算机模拟该随机现象,通过对大量模拟仿真试验的结果来分析计算所求参数,得出实际问题的近似解.4晋中学院JinzhongUniversity晋中学院JinzhongUniversity蒙特卡罗方法的特点可归纳成三个方面:(1)蒙特卡罗方法及其程序结构简单.(2)蒙特卡罗方法的收敛性及收敛速度与问题的维数无关.(3)蒙特卡罗方法的适用性强,可用在求很多解析方法或常规数值方法难解问题的低精度解.5晋中学院JinzhongUniversity晋中学院JinzhongUniversity1.1.1蒲丰投针问题著名的投针问题是几何概率一个早期的例子,它是由法国科学家蒲丰(Buffon)在1777年提出的,因而被称之为蒲丰投针问题.蒲丰投针问题的解决不仅较典型的反映了几何概率的特征及处理方法,而且还可以由此了解蒙特卡洛(Monte一Carlo)方法.6晋中学院JinzhongUniversity晋中学院JinzhongUniversity蒲丰投针问题:平面上画有等距离的平行线,每两条平行线之间的距离为d,向平面任意投掷一枚长为lld的针,试求针与平行线相交的概率.解：设x表示针落下后针的中点M到最近的一条平行线的距离,表示针与平行线所成的角(见图1.1.1),则0,02dx.而针与一直线相交的充要条件是：sin2lx.7晋中学院JinzhongUniversity晋中学院JinzhongUniversity我们把及x表为平面上一点的直角坐标,则所有基本事件可以用边长为及2d的矩形内的点表示出来,而“针与直线相交”这一事件所包含的基本事件可以用图1.1.2中阴影部分内的点表示出来,而所求概率02sin22lldPdd.8晋中学院JinzhongUniversity晋中学院JinzhongUniversity晋中学院JinzhongUniversity晋中学院JinzhongUniversity(1)取白纸一张,在上面画许多间距为d的等距平行线；(2)取一根长度为lld的均匀直针,随机地向画有平行线的纸投去,共投n次,观察针和直线相交的次数m；(3)直线与针相交概率P的近似值可用mn得到,其中n是总投针次数而m是与任一直线相交的总次数,进而由可得的近似值为2nlmd.9晋中学院JinzhongUniversity晋中学院JinzhongUniversity浦丰投针实验模拟程序(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;y9=18;plot(x1,y1,'-b',x1,y2,'-b',x1,y3,'-b',x1,y4,'-b',x1,y5,'-10晋中学院JinzhongUniversity晋中学院JinzhongUniversityb',x1,y6,'-b',x1,y7,'-b',x1,y8,'-b',x1,y9,'-b','LineWidth',2)whilek=nholdonxi=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);11晋中学院JinzhongUniversity晋中学院JinzhongUniversityx=[xi,xj];y=[yi,yj];zz=0.5*l*sin(theta);axis([020020])h_point1=plot(x,y,'-g','LineWidth',1);holdonh_point2=plot(xx,yy,'.r','LineWidth',1);pause(0.01)%暂停0.5秒set(h_point1,'color','k')%将当前直线的颜色由原来的颜色变为黑色k=k+1;12晋中学院JinzhongUniversity晋中学院JinzhongUniversityframe(k)=getframe(gcf);endwritegif('Needlecastexperiment.gif',frame,0.1);所调用的函数程序如下:functionres=writegif(name,frames,dt)nframe=length(frames);forii=1:nframe[image,map]=frame2im(frames(ii));[im,map2]=rgb2ind(image,128);ifii==1IMWRITE(im,map2,name,'GIF',…13晋中学院JinzhongUniversity晋中学院JinzhongUniversity'overwrite','DelayTime',dt,'LoopCount',inf);elseIMWRITE(im,map2,name,…'WriteMode','append',…'DelayTime',dt);%,'LoopCount',inf);endend设2,1dl,(总投针数n=10000000).Matlab实际计算程序(CH1.1.2)：clear14晋中学院JinzhongUniversity晋中学院JinzhongUniversityd=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;15晋中学院JinzhongUniversity晋中学院JinzhongUniversityy=0.5*l*sin(xi);%产生曲线上对应xi的函数值yforii=1:nifyi(ii)=y(ii)m=m+1;%统计处在曲边梯形内的随机点数目endendPi=n/m共计算10次,其结果见表1.1.1.16晋中学院JinzhongUniversity晋中学院JinzhongUniversity表1.1.1计算次数12345计算值3.14433.14193.14273.14413.1426计算次数678910计算值3.14283.13993.14303.14253.14141.1.2计算圆周率的其它方法著名的投针问题利用单位圆与边长为1的正方形面积之比来计算的近似值.具体思想如下:17晋中学院JinzhongUniversity晋中学院JinzhongUniversity如图1.1.4所示,单位圆的14为一个扇形G,它是边长为1的正方形的一部分.考虑扇形面积在正方形面积中所占的比例k,得出其结果为4,然后乘以4就可以得到的值.这里如何计算比例k,运用蒙特卡罗方法的随机投点思想.在正方形中随机投入很多点,使所投点几何概率在正方形中每一个位置的机会均等,然后考察有多少点落在扇形内.其中落在扇形内的点的个数m与18晋中学院JinzhongUniversity晋中学院JinzhongUniversity投点总数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);19晋中学院JinzhongUniversity晋中学院JinzhongUniversityholdonx1=linspace(0,1,500);y1=1;x2=1;y2=linspace(0,1,500);plot(x1,y1,'-k',x2,y2,'-k','LineWidth',2)axis([0101])set(gca,'XTick',0:1)holdon%模拟投掷实验k=0;20晋中学院JinzhongUniversity晋中学院JinzhongUniversitywhilek=nholdonxi=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);21晋中学院JinzhongUniversity晋中学院JinzhongUniversityendwritegif('Throwingexperiment.gif',frame,0.1);%计算pi值xi=unifrnd(0,1,n,1);yi=unifrnd(0,1,n,1);%产生n个符合题意随机点的坐标(xi,yi)m=0;forii=1:nifyi(ii)^2+xi(ii)^2=1m=m+1;%统计处在曲边梯形内的随机点数目22晋中学院JinzhongUniversity晋中学院JinzhongUniversityendendmm=mnn=nPi=m/n*4所调用的函数程序如浦丰投针实验.Matlab实际计算程序(CH1.1.4)：cleard=1;n=10000000;23晋中学院JinzhongUniversity晋中学院JinzhongUniversityxi=unifrnd(0,1,n,1);yi=unifrnd(0,1,n,1);%产生n个符合题意随机点的坐标(xi,yi)k=0;forii=1:nifyi(ii)^2+xi(ii)^2=1k=k+1;%统计处在曲边梯形内的随机点数目endendPi=k/n*4实际计算投掷次数取n=10000000,计算结果见表24晋中学院JinzhongUniversity晋中学院JinzhongUniversity1.1.2.表1.1.2计算次数12345计算值3.14213.14153.14173.14193.1423计算次数678910计算值3.14163.14103.14103.14123.1416从上述结果来看,计算的结果涨落比较大,如果要达到比较高的精度,需要大幅增加投掷的次数.25晋中学院JinzhongUniversity晋中学院JinzhongUniversity1.2利用蒙特卡罗法求解数值积分和函数极值1.2.1定积分的蒙特卡罗算法假定函数fx在,ab内有界连续且0fx,计算定积分baIfxdx.26晋中学院JinzhongUniversity晋中学院JinzhongUniversity为计算出定积分值,可构造一概率模型:取一个边长分别为ba和h的矩形D(如图1.2.1),使曲边梯形在矩形域之内,并在矩形内随机投点,假设随机点均匀地落在整个矩形之内,则落在图中灰色区域内的随机点数k与投掷点总数N之比kN,就近似地等于灰色区面积与矩形面积之比,从而得出定积分kIbahN.27晋中学院JinzhongUniversity晋中