计算物理初步.pdf

《计算物理初步.pdf》由会员分享，可在线阅读，更多相关《计算物理初步.pdf（90页珍藏版）》请在得力文库 - 分享文档赚钱的网站上搜索。

1、计算物理初步晋中学院物理与电子工程学院目录第 1章 蒙特卡罗方法的应用.11.1 蒙特卡罗方法简介.11.1.1 蒲丰投针问题.11.1.2 计算圆周率的其它方法.31.2 利用蒙特卡罗法求解数值积分.61.2.1 定积分的蒙特卡罗算法.61.2.2 二重积分的蒙特卡罗算法.71.2.3 蒙特卡罗求一元函数的最值.71.3 基于蒙特卡罗的电子双缝衍射的计算机模拟.91.3.1 电子双缝衍射的随机模拟.91.3 3 电子双缝衍射图.11第 2 章有限差分法.202.1 有限差分法基础.202.2 拉普拉斯方程.232.2.1 二维场域的边界条件.242.2.2 简单迭代法.242.2.3 超松弛

2、法.252.2.4 应用举例与计算程序.252.3 本征值方程.332.3.1-维薛定四方程的有限差分解法.332.3.2 二维薛定谓方程的有限差分解法.492.4 热传导方程.602.4.1 热传导方程概述.602.4.2 显式差分格式.612.4.3 隐式六点差分格式（CN 格式）.642.4.4 二维扩散方程的有限差分法.722.5 波动方程.802.5.1 显式差分格式.802.5.2 初始、边界条件的差分格式.802.5.3 计算示例.82第1章蒙特卡罗方法的应用1.1蒙特卡罗方法简介蒙特卡罗方法也称随机模拟方法，有时也称作随机抽样技术或统计实验方法.它的基本思想是：首先建 立 个概

3、率模型或随机过程,使它的参数等于问题的解;然后利用计算机模拟该随机现象，通过对大量模拟仿真试验的结果来分析计算所求参数，得出实际问题的近似解.蒙特卡罗方法的特点可归纳成三个方面：(1)蒙特卡罗方法及其程序结构简单.(2)蒙特卡罗方法的收敛性及收敛速度与问题的维数无关.(3)蒙特卡罗方法的适用性强，可用在求很多解析方法或常规数值方法难解问题的低精度解.1.1.1蒲丰投针问题著名的投针问题是几何概率一个早期的例子，它是由法国科学家蒲丰(Buffon)在 1777年提出的,因而被称之为蒲丰投针问题.蒲丰投针问题的解决不仅较典型的反映了几何概率的特征及处理方法,而且还可以由此了解蒙特卡洛(Monte

4、Carlo)方法.蒲丰投针问题：平面上画有等距离的平行线，每两条平行线之间的距离为，向平面任意投掷一枚长为/(/d)的针，试求针与平行线相交的概率.解：设x 表示针落下后针的中点M 到最近的一条平行线的距离，夕表示针与平行线所成的角(见下 图 1.1.1),则0 x y,0 4*4 万 (1.1.1)而针与一直线相交的充要条件是我们把3及 X表为平面上一点的直角坐标，则所有基本事件可以用边长为 及d/2 的矩形内的点表示出来，而“针与直线相交”这一事件所包含的基本事件可以用图1.1.2中阴影部分内的点表示出来，而所求概率I(1.1.3)(1)取白纸一张，在上面画许多间距为 的 等距平行线；(2

5、)取一根长度为/(/d)的均匀直针，随机地向画有平行线的纸投去，共投”次，观察针和直线相交的次数?；浦丰投针实验模拟程序(cxl)：cleard=2;%平行直线间隔1=1;%针长k=0;n=100;%投掷次数xl=O:0.01:20;y 1 =2;y2=4;y3=6;y4=8;y5=10;y6=12;y7=l 4;y8=16;y9=18;plot(x l,y 1 ；-b；x l,y2；-b,x l,y3；-b;x l,y4；-bx l,y5；-b；x l,y6；-b,x l,y7；-b,x l,y8；-bx 1 ,y9；-b7LineWidth,2)while k=nhold onxi=uni

6、fmd(2,18,1,1);yi=unifrnd(2,18,l,l);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 201)h_point 1 =plot(x,y,-g,/LineWidth 1);hold onh_point2=plot(xx,yy,.r,Line Width,1);pause(O.Ol)%暂停0.5 秒set(h_point

7、l；color,；k,)%将当前直线的颜色E ll原来的颜色变为黑色k=k+1;end图1.1.3是投掷100次的结果.(3)直线与针相交概率P的近似值可用m/n得到，其中n是总投针次数而m是与任一直线相交的总次数，进而由可得万的近似值为2 n li md.设d=2,/=1,Matlab实际计算程序如下(总投针数n=10000000)(cx2)：clear2d=2;1=1；theta=pi;n=10000000;xi=unifrnd(0,theta,n,1);yi=unifrnd(0,d/2,n,1);m=0;y=0.5*l*sin(xi);for i=1:nif yi(i)=y(i)m=m+

8、1;end%产牛由个符合题意随机点的坐标(xi,yi)%产生曲线上对应xi的函数值y%统计处在曲边梯形内的随机点数目endPi=n/m共计算六次，其结果见表1.1.1.表1.1.1计算次数12345678910计算值 3.14433.14193.14273.1441 3.14263.14283.13993.14303.14253.14141.1.2计算圆周率的其它方法利用单位圆与边长为1的正方形面积之比来计算的近似值.具体思想如下：如图1.1.4所示，单位圆的1/4为一个扇形G,它是边长为1的正方形的一部分.考虑扇形面积在正方形面积中所占的比例k,得出其结果为;r/4,然后乘以4就可以得到/的

9、值.这里如何计算比例k,运用蒙特卡罗方法的随机投点思想.在正方形中随机投入很多点，使所投点著名的投针问题是几何概率一在正方形中每一个位置的机会均等，然后考察有多少点落在扇形内.其中落在扇形内的点的个数机与投点总数”之 比就是A的近似值.模拟实验程序(cx3)：clear图 1.1.4d=l;n=1000;%绘制几何图形x=linspace(-1,1,1000);y=sqrt(l-x.A2);plot(x,y；r,LineWidth,2);hold onx 1 =linspace(0,1,500);yi=i;3x2=l;y2=linspace(0,1,500);plot(x 1 ,y 1,-kx

10、2,y2,-k,LineWidth*,2)axis(0 1 0 1)set(gca/XTick0:1)hold on%模拟投掷实验k=l;while k=nhold onxi=unifmd(0,1,1,1);yi=unifrnd(0,1,1,1);plot(xi,yi,1)h_point=plot(xi,yi,.r,LineWidth 1);pause(0.0001)%暂停0.0001 秒setCh-point/colof/b1)%将当前点的颜色由原来的颜色变为蓝色k=k+l;end%计算pi值xi=unifrnd(0,l,n,l);yi=unifmd(0,l,nJ);%产生。个符合强意随机点

11、的坐标(xi、yi)m=0;for i=1:nifyi(i)八 2+xi 八 2=1m=m+1;%统订处在曲边梯形内的随机点数目endendmm=mnn=nPi=m/n*4Matlab实际计算程序(cx4)：cleard=l;n=10000000;xi=unifrnd(0,l,n,l);yi=unifrnd(0,l,n,l);%产:生n个符合题意随机点的坐标(xi,yi)k=0;for i=1:nif yi(i)A2+xi(i)A2=lk=k+l;%统计处在曲边梯形内的随机点数目endend4Pi=k/n*4实际计算投掷次数取n=10000000,计算结果见表1.12表1.1.2计算次数123

12、45678910计算值 3.1421 3.14153.14173.14193.14233.14163.14103.14103.14123.141651.2 利用蒙特卡罗法求解数值积分1.2.1 定积分的蒙特卡罗算法假定函数K 0在出,句内有界连续且f(x)o,对定积分/=f/(x)d x,(1.2.1)为计算出定积分值，可构造一概率模型：取一个边长分别为儿。和 h 的矩形Q(如图1 2 1),使曲边梯形在矩形域之内，并在矩形内随机投点，假设随机点均匀地落在整个矩形之内，则落在图中灰色区域内的随机点数k 与投掷点总数N 之比k/N,就近似地等于灰色区面积与矩形面积之比，从而得出定积分I=-(b-

13、a)h (1.2.2)O a b x例 1.2.1:计算 dx.图 1 2 对 x)=e 3,我们很难求出其原函数，所以用牛顿-莱布尼茨公式无法求解(见图1.2.2).这个问题可运用蒙特卡罗方法依如下步骤求出其近似值.步 骤 1:产生矩形域D 内的N 个均匀随机点Pi(xi,yi);步骤2:统计满足条件yiWf(x)的落在曲边梯形内随机点Pi的数目k;步骤3:(取h=l,a=-1,b=1,N=100000),则定积分的近似值I弋(b-a)hk/N=k/5000.Matlab计算程序(cx5)：clear图 1.2.2xi=unifrnd(-1,1,100000,1);yi and(100000

14、,l);%产 生 100000个符合题意随机点的坐标(xi,yi)k=0;y=exp(-xi.A2);%产生曲线上对应x i的函数佗fori=l:100000if yi(i)=y(i)k=k+l;endendI=k/50000%统计处在曲边梯形内的随机点数H将以上程序利用matlab软件重复执行6 次，所得计算结果见表1.2.1.表 1.2.1第 i 次123456积分近似值1.49341.49281.49631.49361.49461.4941注意：山于随机误差的影响，每次计算机模拟的结果可能不一样，但是各次计算结果总是在准确值附近作微小摆动.61.2.2 二重积分的蒙特卡罗算法对二重积分/

15、=,设 x,y)为区域A上 的 有 界 函 数 且2 0 ,据其几何意义,它是以x,y)为曲面顶，A为底的曲顶柱体C的体积.据此，用均匀随机数计算二重积分的蒙特卡罗方法基本思路为：假设曲顶柱体C包含在己知体积为力的几何体。的内部，在。内产生N个均匀随机点，统计出在C内部的随机点数目N,则 Nc.N N例 1.2 2 计算+/_/_/)_,2+.2卜xdy,其中 A =A(x,y +y2 4 1.A利用M a tla b先做出函数的图像(见图1.2.3).图1.2.3分析:该二重积分可看作以z=(1+7 1二7二7)-次17为顶的曲顶柱体的体积，此曲顶柱体在一个边长为2的立方体之内(见图1.2.

16、3),可计算出其精确值为乃.现利用蒙特卡罗算法计算其近似值，见以下程序：M a tla b计算程序(c x6)：c le a rn=4 0 0 0 0 0;x=uni f md(-ly=uni f rnd(-l,l,n,l);zi =uni f md(0,2,n,l);z=l+sq rt(l-x.A2-y.A2)-sq rt(x.A2+y.A2);N c =0;f o r i =1:n图 1.2.3i f x(i)A2+y(i)A2 =1&zi(i)=z(i)N c =N c+l;表 1.2.2e nde ndI =8*N c/n下面是重复十次运算所得数据(见表L 2.2)第k次1234567

17、8910近似值3.14 2 9 3.14 2 63.14 19 3.14 5 33.14 103.14 133.14 2 8 3.13 933.13 923.14 3 8从上面十个数据可以看出，每次计算结果在精确值区的附近摆动.为提高精度，可选取其它随机数促进蒙特卡罗的收敛性，例如选取有利随机数来代替均匀随机数可在较少随机点的情况下使计算精度大大提高.1.2.3 蒙特卡罗求一元函数的最值高数中求一元函数f(x)在闭区间 a,b 的最值的方法为：找出函数f(x)的所有驻点和不可导点并计算出其相应函数值，并将其与区间端点函数值比较得出函数的最大值最小值.但实际问题中方程F(x)=0很难求出其根，所

18、以驻点往往很难求出.若利用蒙特卡罗方法：在 a,b 上产生n个随机数，计7算出这些数的函数值并作比较，则可得出最值近似值.例 1.2.3:f(x)=(l-x3)sin3x,2 五 WxW2 n,求函数的最大值.绘制函数图形(图1.2.4).解法一：使用蒙特卡罗方法Matlab计算程序(cx7)：clearn=le4;x=unifmd(-2*pi,2*pi,n,l);fx=(l-x.A3).*sin(3*x);fmax=max(fx)将上述程序重复运行6 次(见表1.2.3),则可得出函数的最大值为194.9062表 1.2.3第 i 次123456最大值194.906194.9062194.9

19、055194.9056194.9061194.9062解法二：采用等距法搜索函数最值Matlab计算程序(cx8)：clearn=le4;x=linspaceQ2*pi,2*pi,n);fx=(l-x.A3).*sin(3*x);fmax=max(fx)采用等距法搜索函数最值得函数最大值为194.9061.根 据 图 1.2.4,可知两种方法的结果都很接近函数最大值.如果取n=le5所得结果亦为194.9062.81.3基于蒙特卡罗的电子双缝衍射的计算机模拟电子双缝衍射是微观粒子具有波动性的重要证明实验,才由Tonomura等人做出了该实验鉴于一般教学仪器不具备进行上述实验的条件，本文依据电子

20、衍射的几率密度函数，运用蒙特卡罗随机模拟方法，借助计算机的数据可视化技术、绘图技术，构 建 电 子双缝衍射的动态随机过程，清晰地演示出电子衍射的全过程.1.3.1 电子双缝衍射的随机模拟 一电子衍射是大量电子随机运动的统计结果，可开始是作为假想实验而提出的，1988年图1.3 电子双缝衍射示意图运用蒙特卡罗方法对其进行模拟.运用蒙特卡罗方法来处理本文的问题，首先必须根据要处理问题的规律，建立一个概率模型，然后进行随机抽样试验，从而得出一组按已知分布的随机数序列.最后依据这一随机数序列，借助计算机程序设计语言或图形软件，就可实现电子双缝衍射动态随机过程的模拟.1.概率模型的构建图1 3 1是 电

21、 子 双 缝 衍 射 的 示 意 图.衍 射 屏 位 于 平 面，观测屏位于x o y平面，缝$和S2的宽度均为a,两缝中心间距为(a+b),衍射屏到观测屏的距离为D.设动量为p、能量为E的自由电子沿Z轴正方向入射到双缝，其波函数为：(1.3.1)设在t=0时刻，波前到达双缝处(Z=0),由式(1)可知，在双缝处波函数为一常数.为简单起见，设2式中,k=2 w/入.所以，电子经过双缝在观测屏上P点出现的几率密度为:(1.3.2)(1.3.3)sin2 1k sin e)cos2 k2sin (a+b),则有:代入式(1.9)有：(1.3.4)(1.3.5)-1-xADsin。(Ax),.、%-

22、ycos(Bx)(1.3.6),na 乃(a+b)具中 A=x,B=i-AD AD2.数据的生成根据上述电子在观测屏上出现的几率密度函数进行随机抽样，便得到按此几率密度函数分布的随机数序列.在蒙特卡罗方法中，有多种方法可实现按已知分布的随机抽样.根据本文要处理的问题的特点，采用Von Neumann的舍选法.舍选法的具体做法是：(1)计算机在一定的范围内随机地选取观测屏坐标点(冷必)，并计算H=w(xJ/w11m.其中，w(x j是几率密度函数w(x)在点(%)的值，w皿是几率密度函数w(x)的最大值；(2)计算机产生一个0 至 1之间均匀分布的随机数M;(3)将 H与M 进行比较，若”N M

23、,则选取该点(4 力),若H%，所以:叱而(Ar)(1.3.7)根据以上讨论编写的电子双缝衍的Matlab程序如下(cx9):clearh=6.62559e-34;m=9.10908e-31;10e=1.6021e-19;U=1000;D=0.25;rd=h/sqrt(2*m*e*U);A=(pi*a)/(rd*D);B=pi*(a+b)/(rd*D);axis(-5e-5,5e-5,-4e-5,4e-5)set(gca,color,0.122,0.012,0.62)%计算电子波长rd%计算A 和 B%设置坐标轴标度范围%设置坐标面背景颜色title。电子双缝衍射动态随机过程演示Yfontsi

24、ze1 6,color,k)%设置标题i=l；while k=100%用循环结构控制电子的数目x=5e-5*(2*rand-l);%产生5e-5至-5e-5之间均匀分布的随机数赋给xy=4e-5*(2*rand-l);%产生4e-5至-4e-5之间均匀分布的随机数赋给yH=(sin(A*x)A2)/(A*x)A2)*(cos(B*x)A2);%计算 HM=rand;%产 生 0 至 1之间均匀分布的随机数赋给Mif H=M%用分支结构选择符合条件的坐标点hold on%保留节前坐标系中已存在的图形对象h_point=plot(x,y,.r,EraseMode,none,markerSize,1

25、0);%，用乡 颜色显示符合条件的坐标点,并赋句秉值给h_point;i=i+l;pause(O.OOl)%暂停 1 秒set(h_point,color,w)%将当前坐标点山红颜色变为白颜色endend1.3.3电子双缝衍射图运行上述程序便会得到电子双缝衍射动态随机过程的演示，图 1.3.2、图 1.3.3、图 1.3.4分别是电子数为100、1000、2000时的衍射图.从中可看出，当电子数目较少时(图1.3.2),衍射特征不明显,电子表现出的是粒子性；但随着电子数目的增多(如图1.3.3、图 1.3.4),衍射特征逐渐明显，电子表现出了波动性.将衍射结果与根据量子力学原理得出的衍射几率密

26、度分布函数进行比较.可以看出,衍射图中所反映出的电子点的分布，与几率密度分布函数曲线所反映的几率密度大小分布完全一致.图 1.3.5电子数为20 000时的衍射图样与衍射几率密度分布函数曲线对比图(图1.3.4中 markerSize取 2).11M,0电子双缝步时动志随机逋程演示电子双缄忻射动态期机过程演示图 1.3.2 n=1 0 0图 1.3.3 n=1 0 0 0 ;w学 及 改 打;*.t L *.*,.：-M /=Mhold on%产生0 至 1之间均匀分布的随机数赋给M%用分支结构选择符合条件的坐标点%保留当前坐标系中已存在的图形对象h_point=plot(x,y,.r,Era

27、seMode,none,markerSize,10);%J U 勿颜色显示符介条件的坐标点，并赋句秉值给h_point;i=i+l;pause(O.OOl)set(h_point,color,w)endend%暂 停 1秒%将当前坐标点由红颜色变为白颜色图 1.3.3程序：clearh=6.62559e-34;m=9.10908e-31;e=1.6021e-19;U=1000;a=2e-7;b=le-6;D=0.25;rd=h/sqrt(2*m*e*U);A=(pi*a)/(rd*D);%计算电子波长rd13B=pi*(a+b)/(rd*D);%计算 A 和 Baxis(-5e-5,5e-5,

28、-4e-5,4e-5)%设置坐标轴标度范围set(gca；color;0.122,0.012,0.62)%设置坐标面背景颜色title。电子双缝衍射动态随机过程演示;fontsize；16；color,k)%设置标题i=l；while k=1000%用循环结构控制电子的数目x=5e-5*(2*rand-1);%产 生 5e-5至5e-5之间均匀分布的随机数赋给xy=4e-5*(2*rand-l);%产 生 4e-5至-4e-5之间均匀分布的随机数赋给yH=(sin(A*x)A2)/(A*x)A2)*(cos(B*x)A2);%计算 HM=rand;%产生0 至 12 M均匀分布的随机数赋给Mi

29、f H=M%用分支结构选择符介条件的坐标点hold on%保留当前坐标系中已存在的图形对象h_point=plot(x,y,.r,EraseMode,none,markerSize,10);%用红颜色显示符合条件的坐标点,并i=i+l;pause(O.OOl)set(h_point,color,endend图 1.3.4程序：clearh=6.62559e-34;m=9.10908e-31;e=1.6021e-19;U=l000;a=2e-7;b=le-6;D=0.25;rd=h/sqrt(2*m*e*U);A=(pi*a)/(rd*D);B=pi*(a+b)/(rd*D);axis(-5e-

30、5,5e-5,-4e-5,4e-5J)set(gca/color0.122,0.012,0.62)赋句秉值给h_point;%暂 停 1秒%将当前坐标点山红颜色变为白颜色%计算电子波长rd%计算A 和 B%设置坐标轴标度范围%设置坐标面背景颜色title，电子双缝衍射动态随机过程演示？fontsize；16，cok)”k)%设置标题i=l；while=M%用分支结构选择符合条件的坐标点hold on%保留当前坐标系中已存在的图形对象h_point=plot(x,y,.r,EraseMode,none*,markerSize,10);%JIJ 红颜色(u 示符合条件的坐标，点、,并14i=i+l

31、;赋句秉值给h_point;图 1.3.5程序:clearpause(O.OOl)set(h_point,color,w)endend%暂 停 1秒%将当前坐标点由红颜色变为白颜色title，电子双缝衍射动态随机过程演示;fontsize；16；color；k)%设置标题h=6.62559e-34;m=9.10908e-31;e=1.6021e-19;U=1000;a=2e-7;b=1 e-6;D=0.25;rd=h/sqrt(2*m*e*U);A=(pi*a)/(rd*D);%计算电子波长rdB=pi*(a+b)/(rd*D);%计算A 和 Baxis(-5e-5,5e-5,-4e-5,4e

32、-5J)set(gca；color0.122,0.012,0.62)%设置坐标轴标度范围%设置坐标面背景颜色H=(sin(A*x)八 2)/(A*x)八 2)*(cos(B*x)八 2);%计算 H1=1;while i=Mhold on%产生0 至 1之间均匀分布的随机数赋给M%用分支结构选择符合条件的坐标点%保留当前坐标系中已存在的图形对象h_point=plot(x,y,.r,EraseMode,none,markerSize,10);%用红颜色显示符合条件的坐标点,)-赋句秉值给h_point;i=i+l;pause(O.O)set(h_point,color,w1)endend%暂

33、停 1秒%将当前坐标点由红颜色变为白颜色15图 1.1.2程序：clearclearx=linspace(0,pi,500);y=sin(x).*0.5;fin(x,y；g,);hold onx 1 =linspace(0,pi,500);yi=i;x2=pi;y2=linspace(0,l,500);plot(xl,y-k,x2,y2,-k,LineWidlh,2)textCInterpreter/latex1,.String1,*$x=frac l2,sin left(phi right)Position,pi/2.7 0.6,.FontSize16,.,color,1k1)axis(0

34、piO 1)set(gca,XTick,0:pi:pi)set(gca：XTickLaber,O,)set(gca,YTick*,0:1)set(gca；YTickLabel1,0;d/2,)xlabel(fontsize 24 nnphir)ylabel(fontsize 24 rmxr)%end绘制图4的程序：clearx=linspace(-1,1,1000);y=sqrt(l-x.A2);fill(x,y；b,);hold onx l=linspace(0,l,500);yl=l;x2=l;y2=linspace(0,1,500);plot(x l,y l,-k;x2,y2；-k,；L

35、ineWidth,2)textCInterpreter/latex,.,String；$xA2+yA2=l$；.Position,0.5 0.5,.*FontSize,16,.colof/f)axis(0 1 0 1)set(gca；XTick,0:1)16set(gca；XTickLaber,O；r)set(gca,YTick,O:1)set(gca；YTickLabel,1)xlabel(fontsize 24rmx,)ylabel(fontsize 24rmy)%end图 1.6程序:clearfigh=figure(color/w);xl=-l:0.01:1;yl=exp(-xl.A2

36、);x2=-l;y2=0:0.001:1;x3=l;y3=0:0.001:1;x4=-1:0.001:1;y4=i;x5=-1.2:0.001:1.2;y5=0;x6=0;y6=0:0.001:1.2;plot(x l,y l,-k,x2,y2,-k,x3,y3,-k,x4,y4,-k,x5,y5,-k,x6,y6,-k7LineWidth,2)axis(-1.5 1.5 0 1.5)set(gca,*XTick-I:1:1)set(gca；XTickLaber,)setCgca/YTickO:1:1)%end图 1 2 3 程序clearR=l;x 1 =linspace(-1,1,200)

37、;y l=linspace(-l,1,200);x,y=meshgrid(x l,yl);z=l+sqrt(1 -x.A2-y.A2)-sqrt(x.A2+y.A2);r=sqrt(x.A2+y.A2);a,b=find(rR);for i=l:length(a);z(a(i),b(i)=nan;endsurf(x,y,z)shading flat17xlabel(fontsize 14bfx)ylabel(fonlsize 14bfy)zlabel(Vontsizef 14bfz)axis(-l 1 -1 10 2)box on%end图 1.2.4程序：clearx=-2*pi:0.000

38、01:2*pi;y=(1 -x.A3).*sin(3*x);plot(x,y,-k,LineWidth；2);title(fontsize 18bfy=(l-xA3)*sin(3*x)xlabel(Afontsize 16bfx*)ylabel(Afontsizef 16bfy*)axis(-2*pi 2*pi-150 200J)grid on%end图 1.3.5程 序 1：clearh=6.62559e-34;m=9.10908e-31;e=1.6021e-19;U=1000;a=2e-7;b=le-6;D=0.25;rd=h/sqrt(2*m*e*U);A=(pi*a)/(rd*D);B

39、=pi*(a+b)/(rd*D);axis(-5e-5,5e-5,-4e-5,4e-5)set(gca,color,0.122,0.012,0.62)%计算电子波长rd%计算A 和 B%设置坐标轴标度范围%设置坐标面背景颜色title，电子双缝衍射动态随机过程演示？fontsize；16，colo,k)%设置标题i=l；while k=20000%用循环结构控制电子的数目x=5e-5*(2*rand-l);勿产生5e-5至-5e-5之间均匀分布的随机数赋；令xy=4e-5*(2*rand-l);%产牛.*5 至：-4c-5之1 川均匀分布的随机数赋给yH=(sin(A*x)A2)/(A*x)A

40、2)*(cos(B*x)A2);%计算 HM=rand;%心 生 0 至 之间均匀分布的随机数赋给Mif H=M%用分支构选捽行介条件的坐标点hold on%保留当前坐标系中存在的图形对象h_point=plot(x,y,.r,EraseMode,none,markerSize,2);%用红颜色!u 示符介条件的坐标点、并赋句秉值给h_point;18i=i+1;pause(O.O)set(h_point,color1,w)endend图 1.3.5程序2：clearh=6.62559e-34;m=9.10908e-31;e=1.6021e-19;U=1000;a=2e-7;b=1 e-6;D

41、=0.25;rd二 h/sqrt(2*m*e*U);A=(pi*a)/(rd*D);B=pi*(a+b)/(rd*D);x=-5e-4:le-10:5e-4;%暂停0 秒%将当前坐标点由红颜色变为白颜色%计算电子波长rd%计 算 A 和 Bw=(sin(A*x).A2)./(A*x).A2).*(cos(B*x),A2);plot(x,w,-k,LineWidth,0.5)axis(-5e-5,5e-5,0J J)grid onxlabel(fontsize 18rmx)ylabel(Vontsize 18rmw/w_0)%end19第 2 章有限差分法物理学和其他学科领域的许多问题往往在被分

42、析研究之后，都可以归纳为常微分或偏微分方程的求解问题.一般说来，处理一个特定的物理问题，除了需要知道它满足的数学方程外，还应当同时知道这个问题的定解条件，然后才能设计出行之有效的计算方法来求解.有限差分法是一种得到广泛应用的较好的数值求解方法.2.1 有限差分法基础在有限差分方法中，我们放弃了微分方程中独立变量可以取连续值的特征，而关注独立变量离散取值后对应的函数值.但是从原则上说，这种方法仍然可以达到任意满意的计算精度.因为方程的连续数值可以通过减小独立变量离散的取值的间隔，或者通过离散点上的函数值插值计算来近似得到.这种方法是随着计算机的诞生和应用发展起来的.其计算格式和程序的设计都比较直

43、观和简单，因而，它的实际应用已经构成了计算数学和计算物理的重要组成部分.我们在这一章中将介绍这种方法和它的应用举例.有限差分法的具体操作分为两个部分：(1)用差分代替微分方程中的微分，将连续变化的变量离散化，从而得到差分方程组的数学形式；(2)求解差分方程组.在第一步中，我们通过所谓的网络分割法，将函数定义域分成大量相邻而不重合的子区域.通常采用的是规则的分割方式.这样可以便于计算机自动实现和减少计算的复杂性.网络线划分的交点称为节点.若与某个节点尸相邻的节点是定义在场域内的节点，则尸点称为正则节点.若节点产处在定义域外的相邻节点，则P点称为非正则节点.在第二步中，数值求解的关键就是要应用适当

44、的计算方法，求得特定问题在所有这些节点上的离散近似值.设函数为/(x),其独立变量x有一很小的增量小=，则该函数x)的增量为：4/(x)=x +/?)-/(x)(2.1.1)称 为 函 数 的 一 阶 差 分，它与微分不同，因是有限量的差，故称为有限差分.而一阶差分z lf(x)除以增量的商，即为一阶差商：4/(X)=+J 2)Ax h.一阶差分仍是独立变量X的函数，类似地，按(2.1.1)式计算一阶差分，就得到二阶差分的差分,就得到二阶差分片/(力，即4 RX)=4/(X+/Z)-4/(X)(2.1.3)显然，只要上述增量人 很小，差分与微分之间的差异将很小.由于一阶导数也=l i m(2.

45、1.4)dx AX是无限小的微分d=|叱4/(x)除以无限小微分公=?吗 不 商.因 此 如 图2.1.1所示，当用2、320点的量来表示函数/(X)的导数时，可近似地用差分表示为:d /X+/7卜”X)dx Ax即用有限差分zl x)除以增量/x 的商来表示表示函数/(x)的一阶导数时，可近似表达为:hdf(x)dx Ax(向前差分)(2.1.5)也称为差商.同理，当 用 1、2 点的量来dx Ax(向后差分)而 用 1、3 点的量来表示尸(x)时，则有：d “(x)_+dx Ax2h(中心差分)(2.1.7)/(-)(2.1.6)A f +h)公)O x-h x x+h x可见，函数f(x

46、)在 X处的导数可用三种差商近似表示，但哪一种 图2.1.1差分运算原理图精确度更高，则可以通过泰勒公式的展开式来分析.f(x +h)=f(x)+h -+-h2-+-八)dx 2!dx2(2.1.8)h和/=土 +L 2d ：，)+.(2.1.9)V)dx 2!dx2可见，将(21.8)式中的二阶以及二阶以上的项略去后，便可得到竺=/(：“)二/(X),这就dx h是(2.1.5)式.同理，由(2.1.9)式可得到(2.1.6)式，它们都截断于人也。项，而把小项以及更高募次的dx项全部略去.(2.1.7)式相当于把相应的泰勒公式(2.1.8)式及(2.1.9)式联立起来，得到：f(x +h)-

47、f(x-h)=2 h -+-h3 -+-(2.1.10)dx 3!ax当略去二项以及更高幕次的项以后，便得到d 付 x+/(X i)dx 2h即(2.1.7)式.很明显，以上三种差商表达式中，以(2.1.7)式所示的差商的截断误差最小，因为(2.1.7)式略去的是三阶以上的无穷小，而其他的是二阶以上的无穷小.对于二阶导数，同样可用 阶差商的差商来近似表示，即/(x)1 (d _ d 叫 U /f(x +/z)dx2 Ax dx I dx J h、h/(x +/z)一)(x)_/(X+/7)_ 2/3 +/(X-)%I h h)厂(2.1.11)上式相当于把泰勒公式小+/1)+小叫=2个)+/冬

48、+，竽+.截断到二q。项，略去了/项以及更高幕次的项.(2.1.12)偏导数也可以仿照上述方法表示为差商，它用各离散点上函数的差商来近似代替该点的偏导数,将需求解的边值问题转化为一组差分方程，而后根据差分方程组(代数方程组)，解出位于各离散点上的待求函数值，便可得到相应的数值解.21由上可见，有限差分法是以差分原理为基础的一种数值方法，它实质上是将物理上连续域的问题变换为离散系统的问题来求解，也就是用网格状离散化模型上各离散点的数值解来逼近连续场域的真实解的.有限差分法的应用范围很广，不但能求解均匀媒质和不均匀线性媒质中的位场，而且还能求解非线性媒质中的场；它不仅能求解恒定场或似稳场，还能求解

49、时变场.在边值问题的数值方法中，此法是相当简便的，在计算机存储容量容许的情况下，有可能采用较精细的网格，使离散化模型能较精细地逼近真实问题，获得具有足够精度的数值解.222.2 拉普拉斯方程二维拉普拉斯方程可以用有限差分法进行近似计算.首先把求解的区域划分成网格，然后把求解区域内连续的场分布用网格节点上的离散的数值解代替.当然，把网格分得充分细，才能达到足够的精度.网格的划分有不同的方法，我们只讨论正方形的划分方法.如 图 2.2.1 所示，在二维平面场中，我们将区域划分为许多正方形网格，每个网格的边长为(称图2.2 有限差分的正方形网格点为步长)，两组平等线的交点称为网格节点.设0点电位为外

50、,。点网格节点的电位为俗、心、心和%.平现在来推导一个用夕1、死、仍和g表示例的公式.如 图 2.2.1 所示，设正方形每个网格边长为(称为步长)，网格节点外八其上下左右四个节点的电位分别为%.”弘 JT，乐 J,给“.在h充分小的情况下，可以外,为基点进行泰勒展开：0 /+1=R /+-h +-h2+匕dy 2 dy2 6 dy3(p i i =(P i j -h+-h2+77 dy 2 dy2 6 dy3“%+和粤 F汐+.把以上四式相加，得：.1 2(52(p%/+1 +(P i.j-x +(P i-x,j +0 7,i =4(P tjh 获显然，在将以上四式相加的过程中，a的所有奇次方