《计算物理初步-计算机科学与技术.pdf》由会员分享，可在线阅读，更多相关《计算物理初步-计算机科学与技术.pdf（158页珍藏版）》请在taowenge.com淘文阁网|工程机械CAD图纸|机械工程制图|CAD装配图下载|SolidWorks_CaTia_CAD_UG_PROE_设计图分享下载上搜索。
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 章有限差分法.212.1 有限差分法基础.212.2 拉普拉斯方程.242.2.1 二维场域的边界条件.252.2.2 简单迭代法.252.2.3 超松弛法.26
2、2.2.4 应用举例与计算程序.262.3 本征值方程.342.3.1 一维薛定帽方程的有限差分解法.342.3.2 二维薛定将方程的有限差分解法.502.4 热传导方程.612.4.1 热传导方程概述.612.4.2 显式差分格式.622.4.3 隐式六点差分格式(CN 格式).652.4.4 二维扩散方程的有限差分法.732.5 波动方程.812.5.1 显式差分格式.812.5.2 初始、边界条件的差分格式.812.5.3 计算示例.83第 3 章转移矩阵.913.1 转移矩阵方法1.913.2 一维任意势垒的透射系数及几率密度分布.933.3 计算实例.953.4 转移矩阵方法2.10
3、7第 4 章 时域有限差分法.1234.1 时域有限差分法的优点.1234.2 麦克斯韦方程.1234.3 FDTD基本方程.1254.4 良匹配层(APML)吸收边界条件.1314.5 数值稳定性分析.1394.6 FDTD法在计算微结构光纤色散特性中的应用.1404.7时域有限差分法的发展.144第 5 章 频域有限差分法.1475.1 频域有限差分法的基本原理.1475.2 各向异性良匹配层的介绍.1495.3 Compact 2D-FDFD 的应用.1512第 1 章 蒙特卡罗方法的应用1.1蒙特卡罗方法简介蒙特卡罗方法也称随机模拟方法，有时也称作随机抽样技术或统计实验方法.它的基本思
4、想是：首先建立一个概率模型或随机过程,使它的参数等于问题的解:然后利用计算机模拟该随机现象，通过对大量模拟仿真试验的结果来分析计算所求参数，得出实际问题的近似解.蒙特卡罗方法的特点可归纳成三个方面：(1)蒙特卡罗方法及其程序结构简单.(2)蒙特卡罗方法的收敛性及收敛速度与问题的维数无关.(3)蒙特卡罗方法的适用性强，可用在求很多解析方法或常规数值方法难解问题的低精度解.1.1.1蒲丰投针问题著名的投针问题是几何概率一个早期的例子，它是由法国科学家蒲丰(Buffon)在 1777年提出的,因而被称之为蒲丰投针问题.蒲丰投针问题的解决不仅较典型的反映了几何概率的特征及处理方法,而且还可以由此了解蒙
5、特卡洛(Monte-Carlo)方法.蒲丰投针问题：平面上画有等距离的平行线，每两条平行线之间的距离为，向平面任意投掷一枚长为/(/d)的针，试求针与平行线相交的概率.解：设x 表示针落下后针的中点M 到最近的一条平行线的距离，夕表示针与平行线所成的角(见下 图 1.1.1),则0 X -,0(p7l：(1.1.1)而针与一直线相交的充要条件是xw gsin。(1.1.2)我们把8及 x 表为平面上一点的直角坐标，则所有基本事件可以用边长为万及2 的矩形内的点表示出来，而“针与直线相交”这一事件所包含的基本事件可以用图1.1.2中阴影部分内的点表示出来，而所求概率71 2(1)取白纸一张，在上
6、面画许多间距为d的等距平行线；(2)取一根长度为/(/)的均匀直针，随机地向画有平行线的纸投去，共投次，观察针和直线相交的次数机；浦丰投针实验模拟程序(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=l 2;y7=l 4;y8=l 6;y9=l 8;plot(x l,y l,-b*,x 1 ,y2,-b,x l,y3,-b*,x 1 ,y4,-b*,x 1 ,y5,-b,x l,y6/-b,x l,y7,-b,x 1 ,y8；bx l,y9,-b7LineWidth;2
7、)while k=nhold onxi=unifmd(2,18,1J);yi=unifrnd(2,18,1,1);theta=unifmd(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*1 *s in(theta);axis(0 20 0 20)hpointl=plot(x,y,-g7Line Wid th1,1);hold onh_point2=plot(xx,yy,.r7LineWid th 1);pause(O.O
8、l)%暂停 0.5 秒set(h_pointl；color7k,)%将当前直线的颜色由原来的颜色变为黑色k=k+l;end图1.1.3是投掷100次的结果.(3)直线与针相交概率P的近似值可用m/n得到，其中n是总投针次数而m是与任一直线相交的总次数，进而由可得万的近似值为设d =2,/=1,Matlab实际计算程序如下(总投针数n=10000000)(cx2)：clear2d=2;1=1;theta=pi;n=l0000000;xi=unifrnd(0,theta,n,1);yi=unifrnd(0,d/2,n,l);m=0;y=0.5*l*sin(xi);for i=1:nif yi(i)
9、=y(i)m=m+l;endendPi=n/m%产生n个符合题意随机点的坐标(xi,yi)%产生曲线上对应xi的函数值y%统计处在曲边梯形内的随机点数目共计算六次，其结果见表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的正方形的一部分.考虑扇形面积在正方形面积中所占的比例上得出其结果为勿4,然后乘以4
10、就可以得到万的值.这里如何计算比例左，运用蒙特卡罗方法的随机投点思想.在正方形中随机投入很多点，使所投点几何概率在正方形中每一个位置的机会均等，然后考察有多少点落在扇形内.其中落在扇形内的点的个数m与投点总数之比就是k的近似值.模拟实验程序(cx3)：cleard=l;n=1000;x图 1.1.4%绘制几何图形x=linspace(0,l,l 000);y=sqrt(l-x.A2);plot(x,y；iVLineWid th；2);hold onx 1 =linspace(0,1,500);yi=i;x2=l;3y2=linspace(0,1,500);plot(x 1 ,y l；-k,x2
11、,y2；-k；LineWid th,2)axis(0 10 1)set(gca；XTick,0:1)hold on%模拟投掷实验k=l;while k=nhold onxi=unifrnd(0,1,1,1);yi=unifmd(0,1,1,1);plot(xi,yi,.r*)h_point=plot(xi,yi；.r7Line Wid th,1);endpause(O.OOOl)set(h_point,color?b)k=k+l;%暂停0.0001秒%将当前点的颜色由原来的颜色变为蓝色%计算pi值xi=unifmd(0,l,n,l);yi=unifmd(0,l,n,l);%产生n个符合题意随机
12、点的坐标(xi,yi)m=0;for ii=1:nifyi(ii)A2+xi(ii)A2=lm=m+l;%统计处在曲边梯形内的随机点数目endendmm=mnn=nPi=m/n*4Matlab实际计算程序(cx4)：cleard=l;n=10000000;xi=unifrnd(O,l,n,l);yi=unifrnd(O,l,n,l);%产生n个符合题意随机点的坐标(xi,yi)k=0;for ii=1:nifyi(ii)A2+xi(ii)A2 0,计算定积分/=公.(1.2.1)为计算出定积分值，可构造一概率模型：取一个边长分别为儿。和 h 的矩形Q(如图1 2 1),使曲边梯形在矩形域之内，
13、并在矩形内随机投点，假 设 随 机 点 均 匀 地 y”落在整个矩形之内，则落在图中灰色区域内的随机点数k 与投掷点总数N 之比k/N,就近似地等于灰色区面积与矩形面积之比，从而得出定积分I=h-a)h (1.2.2)O a b例 1 2 1:计算/=7 d x.图 12对 x)=e 4,我们很难求出其原函数，所以用牛顿-莱布尼茨公式无法求解(见图1.2.2).这个问题可运用蒙特卡罗方法依如下步骤求出其近似值.步 骤 1:产生矩形域D 内的N 个均匀随机点Pi(xi,yi);步骤2:统计满足条件yiWf(x)的落在曲边梯形内随机点Pi的数目k;步骤3:(取h=l,a=-1,b=1,N=1000
14、00),则定积分的近似值I-a)h k/N=W5000.Matlab计算程序(cx5)：clearxi=unifrnd(-1,1,100000,1);yi=rand(100000,1);%产 生 100000个符合题意随机点的坐标(xi,yi)k=0;y=exp(-xi.A2);%产生曲线上对应x i的函数值yfori=l:100000if yi(i)0,据其几何意义,它是以/(x,y)为曲面顶，/为底的曲顶柱体C的体积.据此，用均匀随机数计算二重积分的蒙特卡罗方法基本思路为：假设曲顶柱体C包含在己知体积为匕 的几何体。的内部，在O内产生N个均匀随机点，统计出在C内部的随机点数目七，则/=幺
15、旌.NN例 1.2.2:计算+一y2)_+y2 d xd y,其中力=/(x,y)H+y2 Wl.利用M atlab先做出函数的图像(见图1.2.3).图1.2.3分析:该二重积分可看作以z=(1+J1-i-4)-J x2+/为顶的曲顶柱体的体积,此曲顶柱体在一个边长为2的立方体之内(见图1.2.3),可计算出其精确值为力.现利用蒙特卡罗算法计算其近似值，见以下程序：M atlab计算程序(cx6)：clearn=400000;x=unifrnd(-l,l,n,l);y=unifrnd(-l,l,n,l);zi=unifrnd(0,2,n,l);z=1 +sqrt(1 -x.A2-y.A2)-
16、sqrt(x.A2+y.A2);Nc=0;fb ri=l:nif x(i)A2+y(i)A2=l&zi(i)=z(i)Nc=Nc+l;endendI=8*Nc/n下面是重复十次运算所得数据(见表1.2.2)表 1.2.2第k次12345678910近似值3.1429 3.1426 3.1419 3.1453 3.1410 3.1413 3.1428 3.1393 3.1392 3.1438从上面十个数据可以看出，每次计算结果在精确值n的附近摆动.为提高精度，可选取其它随机数促进蒙特卡罗的收敛性，例如选取有利随机数来代替均匀随机数可在较少随机点的情况下使计算精度大大提高.1.2.3蒙特卡罗求一元
17、函数的最值高等数学中求一元函数/(%)在 闭 区 间 的 最 值 的 方 法 为：找出函数/(%)的所有驻点和不可导点并计算出其相应函数值，并将其与区间端点函数值比较得出函数的最大值最小值.但实际问题7中方程/(x)=0 很难求出其根，所以驻点往往很难求出.若利用蒙特卡罗方法：在 此“上产生个随机数，计算出这些数的函数值并作比较，则可得出最值近似值.例 1.2.3:/(x)=(1 -x3)sin3x,-2 x 2 ,求函数的最大值.绘制函数图形(图124).解法一：使用蒙特卡罗方法Matlab计算程序(cx7)：clearn=le4;x=unifrnd(-2*pi,2*pi,n,l);fx=(
18、l-x.A3).*sin(3*x);finax=max(fx)将上述程序重复运行6 次(见表1.2.3),则可得出函数的最大值为194.9062表 1.2.3第 i 次123456最大值194.906194.9062194.9055194.9056194.9061194.9062解法二：采用等距法搜索函数法求最值Matlab计算程序(cx8)：clearn=le4;x=linspace(-2*pi,2*pi,n);fx=(l-x.A3).*sin(3*x);finax=max(fx)采用等距法搜索函数最值得函数最大值为194.9061.根 据 图 1.2.4,可知两种方法的结果都很接近函数最大
19、值.如果取n=le5所得结果亦为194.9062.图 1.2.481.3基于蒙特卡罗的电子双缝衍射的计算机模拟电子双缝衍射是微观粒子具有波动性的重要证明实验，开始是作为假想实验而提出的，1988年才由Tonomura等人做出了该实验.鉴于一般教学仪器不具备进行上述实验的条件，本文依据电子衍射的几率密度函数，运用蒙特卡罗随机模拟方法，借助计算机的数据可视化技术、绘图技术，构建电子双缝衍射的动态随机过程，清晰地演示出电子衍射的全过程.1.3.1 电子双缝衍射的随机模拟电子衍射是大量电子随机运动的统计结果，可运用蒙特卡罗方法对其进行模拟.运用蒙特卡罗方图1.3.1电子双缝衍射示意图法来处理本文的问题
20、，首先必须根据要处理问题的规律，建立一个概率模型，然后进行随机抽样试验，从而得出一组按已知分布的随机数序列.最后依据这一随机数序列，借助计算机程序设计语言或图形软件，就可实现电子双缝衍射动态随机过程的模拟.1.概率模型的构建图 1.3.1是电子双缝衍射的示意图.衍射屏位于x。平面，观测屏位于必、平面，缝工和S2的宽度均为“，两缝中心间距为(“+b),衍射屏到观测屏的距离为。.设动量为、能量为E 的自由电子沿z 轴正方向入射到双缝，其波函数为：尹二匕上伊甸(1.3.1)设在/=0 时刻，波前到达双缝处(z=0),由式(1)可知，在双缝处波函数为一常数.为简单起见,设在衍射屏上的波函数为：1z(/
21、,0)=0b+-2-a +xz -x aI 2)2 2f 一 a+勺;+勺;-2 (Q +6),则有:代入式(1.9)有：sin=D(1.3.5)w(x)=|/(x,Z)|2=w0-丁(万。=%sm(4 r)cos2(8x)(1.3.6)(小)-A D )2.数据的生成根据上述电子在观测屏上出现的几率密度函数进行随机抽样，便得到按此几率密度函数分布的随机数序列.在蒙特卡罗方法中，有多种方法可实现按已知分布的随机抽样.根据本文要处理的问题的特点，采用Von Neumann的舍选法.舍选法的具体做法是：(1)计算机在一定的范围内随机地选取观测屏一坐标点(占,%),并计算=七)/叱网.其中，w(x,
22、)是几率密度函数w(x)在点(占,%)的值，/是 几 率 密 度 函 数 w(x)的最大值；(2)计算机产生一个。至 1之间均匀分布的随机数M;(3)将“与 进 行 比 较，若H2M,则选取该点(匕,乂)，若则舍去该点，并重复(1)至(3),重新选择坐标点.通过以上方法，便可以得到一个按w(x)分布的随机数序列.这个随机数序列便是我们实现电子衍射随机运动过程的依据.3.电子双缝衍射动态随机运动过程的实现利用上述生成的按w(x)分布的随机数序列，借助计算机程序设计语言或图形软件，就可以实现电子双缝衍射动态随机过程的模拟.本文采用编写Matlab程序的方法来实现.在编写Matlab程序时,用 Ma
23、tlab的 rand 函数来产生均匀分布的随机数，用 plot函数来实现电子衍射点的显示.1.3.2电子双缝衍射的Matlab程序力取 0=2x10-7?,6=1x10?,。=0.25加，加速电子的电压为U=1000.根据4=、：=可算v(x)=s i n )c o s 2(f ix)(s 7)Wmax Ax)根据以上讨论编写的电子双缝衍的Matlab程序如下(cx9):clearh=6.62559e-34;m=9.10908e-31;10e=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=
24、pi*(a+b)/(rd*D);axis(-5e-5,5e-5,-4e-5,4e-5)set(gca,color1,0.122,0.012,0.62)%计算电子波长rd%计算A 和 B%设置坐标轴标度范围%设置坐标面背景颜色title。电子双缝衍射动态随机过程演示Jfontsizd,16Jcolor：k）%设置标题i=l；while i=M%产生0 至 1之间均匀分布的随机数赋给M%用分支结构选择符合条件的坐标点hold on%保留当前坐标系中已存在的图形对象h_point=plot(x,y,.r7EraseMod e,/none,markerSize,10);%用红颜色显示符合条件的坐标点,
25、并赋句秉值给h_point;i=i+l;pause(O.OOl)set(h_point,color?w)end%暂 停 1秒%将当前坐标点由红颜色变为白颜色end1.3.3电子双缝衍射图运行上述程序便会得到电子双缝衍射动态随机过程的演示，图 1.3.2、图 1.3.3、图 1.3.4分别是电子数为100、1000、2000时的衍射图.从中可看出，当电子数目较少时（图 132）,衍射特征不明显,电子表现出的是粒子性；但随着电子数目的增多（如 图 133、图 134）,衍射特征逐渐明显，电子表现出了波动性.将衍射结果与根据量子力学原理得出的衍射几率密度分布函数进行比较.可以看出,衍射图中所反映出的
26、电子点的分布，与几率密度分布函数曲线所反映的几率密度大小分布完全一致.图 1.3.5电子数为20 000时的衍射图样与衍射凡率密度分布函数曲线对比图（图 1 3 4 中 markerSize取 2）.IIx 1n-电子双缝衍射动态随机过程演示5 4 G-2-1 0 1 2 3 4 5图 1.3.5衍射图样与衍射几率密度分布函数曲线对比图12附作图程序：图 1.1.2程序：clearclearx=linspace(0,pi,500);y=sin(x).*0.5;fill(x,y；kf);hold onx 1 =linspace(0,pi,500);yi=i;x2=pi;y2=linspace(0
27、,1,500);plot(xl,yl,-k,x2,y2,-k,LineWid th,2)text(,Interpreter,latex.String7$x=frac 12,sin left(phi right)$.Position1,pi/2.7 0.6,.FontSize16,.lolo”k)axis(0 pi 0 1)set(gca,XTick,0:pi:pi)set(gca,XTickLabel,O?兀 )set(gca,YTick0:1)set(gca；YTickLaber,l0,；d/2,)xlabeKfbntsize 24 rmphT)ylabel(1fbntsize 24 rmx
28、f)%end图1.1.4的程序：clearx=linspace(-1,1,1000);y=sqrt(l-x.A2);fill(x,y；b);hold onx 1 =linspace(0?1,500);yi=i;x2=l;y2=linspace(0,1,500);plot(x l,y 1 ,-kx2,y2,-k7LineWid th,2)text(,Inteq?reter,latex,.String?$xA2+y八 2=1Position1,0.5 0.5,.FontSize,16,.lolorT)axis(0 1 0 1)13set(gca；XTick;O:1)set(gca；XTickLab
29、el;,O,；1*)set(gca,YTick,O:1)set(gca；YTickLabel;,077)xlabelffbntsize 24 rmx*)ylabel(,fbntsize24rmy,)%end图 1.2.2程序：clearfigh=figure(,color,w,);xl=-l:0.01:1;yl=exp(-xl.A2);x2=-l;y2=0:0.001:1;x3=l;y3=0:0.001:1;x4=-l:0.001:1;y4=l;x5=-1.2:0.001:1.2;y5=o；x6=0;y6=0:0.001:1.2;plot(x 1 ,y 1 ,-k,x2,y2,-kx3,y3,
30、-kx4,y4,-kx5,y5,-k,x6,y6,-k7LineWid th,2)axis(-1.5 1.5 0 1,5)set(gca/XTick,-l:1:1)set(gca；XTickLabel,f-1 /O；V)set(gca/YTickO:1:1)%end图 1 2 3 程序clearR=l;x 1 =linspace(-l,l,200);y 1 =linspace(-1,1,200);x,y=meshgrid(x l,yl);z=1 +sqrt(1 -x.A2-y.A2)-sqrt(x.A2+y.A2);r=sqrt(x.A2+y.A2);a,b=find(rR);for i=l:
31、length(a);z(a(i),b(i)=nan;endsurf(x,y,z)shad ing flatxlabel(1fbntsize 14bfx)ylabel Cfontsize 14bfy*)zlabel(fontsize14bf/)axis(-l 1 -1 1 02)box on%end图 1.2.4程序：clearx=-2*pi:0.00001:2*pi;y=(l-x.A3).*sin(3*x);plot(x,y,-k7LineWid th2);titleffontsize 18bfy=(l-xA3)*sin(3*xy)xlabel fd ntsize16bfxf)ylabel(1
32、fbntsize 16bfy)axis(-2*pi2*pi-150 200)grid on%end图 1.3.2程序：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=pi*(a+b)/(rd*D);axis(-5e-5,5e-5,-4e-5,4e-5)setCgca/colofJO.l 22,0.012,0.62)titleC电子双缝衍射动态随机过程演示;fbntsize，,16,color，#)i=l；while i=M
33、%用分支结构选择符合条件的坐标点hold on%保留当前坐标系中已存在的图形对象h_point=plot(x,y,.r,EraseMod e,none7niarkerSize,10);%用红颜色显示符合条件的坐标点，并15赋句秉值给h_point;i=i+l;pause(O.OOl)set(h_point,color?w)end%暂 停 1秒%将当前坐标点由红颜色变为白颜色end图 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*
34、a)/(rd*D);B=pi*(a+b)/(rd*D);axis(-5e-5,5e-5,-4e-5,4e-5)setCgca/colofO.l 22,0.012,0.62)%计算电子波长rd%计算A 和 B%设置坐标轴标度范围%设置坐标面背景颜色title。电子双缝衍射动态随机过程演示Tfbntsize16,color；宴）%设置标题i=l；while i=Mhold on%产生0 至 1之间均匀分布的随机数赋给M%用分支结构选择符合条件的坐标点%保留当前坐标系中已存在的图形对象h_point=plot(x,y,.r,EraseMod e,none1,markerSize,10);%用红颜色显
35、示符合条件的坐标点，并赋句秉值给h_point;i=i+l;pause(O.OOl)set(h_point,color?w)endend%暂 停 1 秒%将当前坐标点由红颜色变为白颜色图 1.3.4程序：clearh=6.62559e-34;m=9.10908e-31;e=1.6021e-19;U=1000;16a=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-5,5e-5,-4e-5,4e-5)set(gca,color1,0.122,0.012,0.62)%计算电子波长rd%
36、计算A 和 B%设置坐标轴标度范围%设置坐标面背景颜色title。电子双缝衍射动态随机过程演示？fbntsize，,color?/)%设置标题i=l；while i=Mhold on%产生0 至 1之间均匀分布的随机数赋给M%用分支结构选择符合条件的坐标点%保留当前坐标系中已存在的图形对象h_point=plot(x,y,.r,EraseMod e,none7markerSize,10);%用红颜色显示符合条件的坐标点,并赋句秉值给h point;i=i+l;pause(O.OOl)set(h_poinLcolorw)endend%暂 停 1秒%将当前坐标点由红颜色变为白颜色图 1.3.5程序
37、：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=pi*(a+b)/(rd*D);axis(-5e-5,5e-5,-4e-5,4e-5)set(gca,color1,0,122,0.012,0.62)%计算电子波长rd%计 算 A 和 B%设置坐标轴标度范围%设置坐标面背景颜色title。电子双缝衍射动态随机过程演示？fontsize；16；color；k)%设置标题i=l；while i=Mhold on%产生0 至 1
38、之间均匀分布的随机数赋给M%用分支结构选择符合条件的坐标点%保留当前坐标系中已存在的图形对象h_point=plot(x,y,.r,EraseMod e,none,markerSize,10);%用红颜色显示符合条件的坐标点，并赋句秉值给h_point;i=i+l;pause(O.O)set(h_point,color?w)endend%暂 停 1秒%将当前坐标点由红颜色变为白颜色图 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=
39、(pi*a)/(rd*D);B=pi*(a+b)/(rd*D);axis(-5e-5,5e-5,-4e-5,4e-5)set(gca/color0,122,0.012,0.62)%计算电子波长rd%计算A 和 B%设置坐标轴标度范围%设置坐标面背景颜色title(电子双缝衍射动态随机过程演示？fontsize，,16 color?k)%设置标题i=l；while i=Mhold on%产生0 至 1之间均匀分布的随机数赋给M%用分支结构选择符合条件的坐标点%保留当前坐标系中已存在的图形对象h_point=plot(x,y，r?EraseMod e?none?markerSize；2);%用红颜
40、色显示符合条件的坐标点，并赋句秉值给h_point;i=i+l；pause(O.O)set(h_point,color,endend%暂停0 秒%将当前坐标点由红颜色变为白颜色图 1.3.5程序2：clear18h=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);%计算电子波长 rdA=(pi*a)/(rd*D);B=pi*(a+b)/(rd*D);%计算 A 和 Bx=-5e-4:le-10:5e-4;w=(sin(A*x).A2)./(A*x).A2).*(cos(B
41、*x).A2);plot(x,w,k?LineWid th,0.5)axis(-5e-5,5e-5,0,l)grid onxlabel(fontsize 18rmx)ylabel(fbntsize 18rmw/w_0)%end1920第2章有限差分法物理学和其他学科领域的许多问题往往在被分析研究之后，都可以归纳为常微分或偏微分方程的求解问题.一般说来，处理一个特定的物理问题，除了需要知道它满足的数学方程外，还应当同时知道这个问题的定解条件，然后才能设计出行之有效的计算方法来求解.有限差分法是一种得到广泛应用的较好的数值求解方法.2.1有限差分法基础在有限差分方法中，我们放弃了微分方程中独立变量
42、可以取连续值的特征，而关注独立变量离散取值后对应的函数值.但是从原则上说，这种方法仍然可以达到任意满意的计算精度.因为方程的连续数值可以通过减小独立变量离散的取值的间隔，或者通过离散点上的函数值插值计算来近似得到.这种方法是随着计算机的诞生和应用发展起来的.其计算格式和程序的设计都比较直观和简单，因而，它的实际应用已经构成了计算数学和计算物理的重要组成部分.我们在这一章中将介绍这种方法和它的一些应用.有限差分法的具体操作分为两个部分：(1)用差分代替微分方程中的微分，将连续变化的变量离散化，从而得到差分方程组的数学形式；(2)求解差分方程组.在第一步中，我们通过所谓的网络分割法，将函数定义域分
43、成大量相邻而不重合的子区域.通常采用的是规则的分割方式.这样可以便于计算机自动实现和减少计算的复杂性.网络线划分的交点称为节点.若与某个节点尸相邻的节点是定义在场域内的节点，则尸点称为正则节点.若节点P处在定义域外的相邻节点，则P点称为非正则节点.在第二步中，数值求解的关键就是要应用适当的计算方法，求得特定问题在所有这些节点上的离散近似值.设函数为/(x),其独立变量x 有一很小的增量小=人，则该函数/(x)的增量为：(x)=x +/O-f(x)(2.1.1)称为函数/(x)的一阶差分，它与微分不同，因是有限量的差，故称为有限差分.而一阶差分4/、(x)除以增量人的商，即为一阶差商：(x)=/
44、(x +)/(x)2)Ax h一阶差分仍是独立变量x的函数，类似地，按(2.1.1)式计算一阶差分，就得到二阶差分的差分,就得到二阶差分d/(x),即=J/(X +/?)-4/(X)显然，只要上述增量 很小，差分/(X)与微分之间的差异将很小.由于一阶导数也%i m 2 3dx 4rT。AX是无限小的微分力。)=l i m 4 f(x)除以无限小微分公=l i m/x 商.因此 如图2.1.1所示，当用2、3(2.1.3)(2.L4)21点的量来表示函数/(X)的导数时，可近似地用差分表示为:小L 2 3=/(x+)-x)(向前差分)dx Ax h(2.1.5)即用有限差分4 f(x)除 以
45、增 量 的 商 来 表 示 也 ，23 也称为差商.同理，当用1、2点的量来dx Ax表示函数/(X)的一阶导数时，可近似表达为：力(X)J/(x)/(x)-/(x-A)-=-=-dx Ax而 用 1、3点的量来表示/(x)时，则有：df(-V)4 x)/(x +A)-/(x-A)-J S S -dx Ax(向后差分)(2.1.6)(中心差分)(2.1.7)图2.1.1差分运算原理图2/7可见，函数/(X)在 X处的导数可用三种差商近似表示，但哪一种精确度更高，则可以通过泰勒公式的展开式来分析.x +)=/(x)+心 詈+仪学+(2.1.8)和/.(X i)=/(x)-力 竽+】小乌国+(2.
46、1.9)ax 21 ax可见，将(2.1.8)式中的二阶以及二阶以上的项略去后，便可得至I J 必1+这就dx h是(2.1.5)式.同理，由(2.1.9)式可得到(2.1.6)式,它们都截断于也且项，而把/项以及更高嘉次的dx项全部略去.(2.1.7)式相当于把相应的泰勒公式(2.1.8)式及(2.1.9)式联立起来，得到：/(x +Q-j)=2/7 +1；小I+(2.1.10)当略去/P 项以及更高幕次的项以后，便得到df(x)+dx 2h即(2.1.7)式.很明显，以上三种差商表达式中，以(2.1.7)式所示的差商的截断误差最小，因为(2.1.7)式略去的是三阶以上的无穷小，而其他的是二
47、阶以上的无穷小.对于二阶导数，同样可用一阶差商的差商来近似表示，即d2f(x)1 f#(x)|#(x)|ifA f(x +h)dx2 dx L dx J 力 1 h1 ”(x +A)-/(x)/(x +/?)-2/(x)+/(%-/?)-=-(2.1.11)上式相当于把泰勒公式/(x +)+/(x-/7)=2/(x)+2+1 4 i l +i (2.1.12)“人 I .”人截 断 到 后 丝 立 项，略去了川项以及更高嘉次的项.dx偏导数也可以仿照上述方法表示为差商，它用各离散点上函数的差商来近似代替该点的偏导数,将需求解的边值问题转化为一组差分方程，而后根据差分方程组(代数方程组)，解出位
48、于各离散点上22的待求函数值，便可得到相应的数值解.由上可见，有限差分法是以差分原理为基础的一种数值方法，它实质上是将物理上连续域的问题变换为离散系统的问题来求解，也就是用网格状离散化模型上各离散点的数值解来逼近连续场域的真实解的.有限差分法的应用范围很广，不但能求解均匀媒质和不均匀线性媒质中的位场，而且还能求解非线性媒质中的场；它不仅能求解恒定场或似稳场，还能求解时变场.在边值问题的数值方法中，此法是相当简便的，在计算机存储容量容许的情况下，有可能采用较精细的网格，使离散化模型能较精细地逼近真实问题，获得具有足够精度的数值解.232.2拉普拉斯方程二维拉普拉斯方程可以用有限差分法进行近似计算
49、.首先把求解的区域划分成网格，然后把求解区域内连续的场分布用网格节点上的离散的数值解代替.当然，把网格分得充分细，才能达到足够的精度.网格的划分有不同的方法，我们只讨论正方形的划分方法.如 图221所示，在二维平面场中，我们将区域划分为许多正方形网格，每个网格的边长为（称6图2 2 1有限差分的正方形网格点为步长），两组平行线的交点称为网格节点.设。点电位为弘，。点网格节点的电位为百、夕2、仍和.平现在来推导一个用夕I、0、0和%表示9 0的公式.如 图2.2.1所示，设正方形每个网格边长为人（称为步长），网格节点外八其上下左右四个节点的电位分别为例仆”/JT，弘 一 叮.在h充分小的情况下，
50、可以”,为基点进行泰勒展开：0川=%0-1 J%T =%M 1 3 V/2+h+-三 h+%3+dy 2 dy26+一弛h+理/2 dx2二马6加+燃h+L雪2+外+dx 2 d x 6加把以上四式相加，得：势J+I +*J T +9一|J +（PM,）=4 g J +h2（需 +黑 （2.2.1）显然，在将以上四式相加的过程中，的所有奇次方项都抵消了.所以，（4.1 8）式是在略去/及其以上项所得到的，其精度为的二次项.由于场中任何点亿）都满足泊松方程*=答+*=F（兑刃d x d y-式中，爪x,y）为场源，则（2.2.1）式变为：对于无源场，尸（x,y）=0,则二维拉普拉斯方程的有限差分