关闭

关闭

关闭

封号提示

内容

首页 计算物理讲义.pdf

计算物理讲义.pdf

计算物理讲义.pdf

上传者: amiao 2010-06-10 评分 0 0 0 0 0 0 暂无简介 简介 举报

简介:本文档为《计算物理讲义pdf》,可适用于自然科学领域,主题内容包含目录i目录第一章基本数学运算第一节插值第二节拟合第三节数值微分第四节数值积分第五节求根第二章常微分方程的初始问题第一节几种简单的数值解法第二节Run符等。

目录i目录第一章基本数学运算第一节插值第二节拟合第三节数值微分第四节数值积分第五节求根第二章常微分方程的初始问题第一节几种简单的数值解法第二节RungeKutta方法第三节多步法第四节稳定性第五节动力学中的有序和混沌第三章边值问题和本征值问题第一节物理学中出现的边值问题和本征值问题举例第二节Numerov算法第三节边值问题的Green函数法第四节打靶法第五节一维Schrödinger方程的定态解第四章特殊函数和Gauss求积法第一节特殊函数第二节Gauss求积法第三节量子散射的Born近似和程函近似第五章矩阵中的数值计算方法第一节物理学中的矩阵第二节矩阵的基本运算第三节一般矩阵的本征值问题第四节对称矩阵的本征值问题第六章椭圆型微分方程第一节离散化和变分原理第二节求解边值问题的一种迭代方法第三节关于离散化的进一步讨论第四节二维定态流体力学第七章抛物型微分方程第一节简单的离散化和条件稳定性第二节隐式格式和回代法第三节高维扩散第四节本征值问题的迭代方法第五节含时间的Schrödinger方程目录ii第八章MonteCarlo方法第一节MonteCarlo方法的基本思想第二节具有特定分布的随机变量得产生第三节Metropolis等人的算法第一章基本数学运算本章将介绍数值计算中一些基本的运算例如插值、拟合、数值微分、数值积分和求根等.第一节插值(Interpolation)当我们要从一组不完全的或离散的数据中获取某些局部的信息时需要使用函数的插值.也就是说对于函数()yfx=(通常是一个未知的或是一个较复杂的函数)xab,若已知ab,上一系列点xxxxn,,,,L处的函数值()yfxii=in=,,,,L当我们要求这些点xxxxn,,,,L之间的点x处的函数值时就需用一个较简单的且满足条件()ϕxyii=in=,,,,L的函数()ϕx来代替()fx这一方法被称为插值法其中()ϕx被称为插值函数()fx被称为被插值函数点ix被称为插值点.当这一较简单的函数()ϕx被选择为多项式时则称为代数插值.当这一较简单的函数()ϕx被选择为三角函数时则称为三角插值.下面主要来介绍代数插值.()线性插值(也称一次插值)已知函数表xxxyyy去构造一个线性多项式()ϕx使其满足()ϕxyii=i=,.这样的一次插值多项式具有如形式()()()ϕxylxylx=其中()lxxxxx=()lxxxxx=被称为线性插值的基函数.误差为()()()()RxfxxOh==ϕ其中xxh=为步长.()二次插值已知函数表xxxxyyyy去构造一个二次多项式()ϕx使其满足()ϕxyii=i=,,.这样的二次插值多项式为()()()()xlyxlyxlyx=ϕ其中()()()()()lxxxxxxxxx=()()()()()lxxxxxxxxx=()()()()()lxxxxxxxxx=被称为二次插值的基函数.误差为()()()()RxfxxOh==ϕ其中h为步长.()n次插值已知函数表xxxxxnyyyyyn去构造一个n次多项式()ϕnx使其满足()ϕniixy=in=,,,,L.n次插值多项式为()()ϕniniixylx==其中()()()()()()()()()()()lxxxxxxxxxxxxxxxxxxxxxiiiniiiiiiin=LLLLin=,,,,L被称为n次插值的基函数.误差为()()()()RxfxxOhnn==ϕ其中h为步长.下面是为n次插值的基函数而编写的MatLat函数文件(chapinterpolybm)functionyy=interpolyb(x,x,n)yy()=fori=:nyy()=yy()*(xx(i))(x()x(i))endyy(n)=fori=:nyy(n)=yy(n)*(xx(i))(x(n)x(i))endforj=:nyy(j)=fori=:jyy(j)=yy(j)*(xx(i))(x(j)x(i))endfori=j:nyy(j)=yy(j)*(xx(i))(x(j)x(i))endend和n次插值函数的函数文件(chapinterpolypm)functiony=interpolyp(x,y,x,n)a=interpolyb(x,x,n)y=a*y'说明当=n和=n时分别可以得到一次插值函数和二次插值函数.例在某实验中要测某短时光束的发光强度随时间的变化已知在实验中测得的一组数据如下t为测数据的时间y为反映光束强度的指标求=t秒处的光的强度指标.t=y=用六次插值多项式来求解该问题.运行的程序为(chapfigurem)x=y=n=min(length(x),length(y))n=nm=h=(max(x)min(x))mfori=:mx(i)=min(x)(i)*hphin(i)=interpolyp(x,y,x(i),n)endninterpolyp(x,y,,n)plot(x,phin)holdonplot(x,y,'bo')该程序的运行结果不仅可以显示出=t秒处的光强度为同时可以画出光强度随时间的演化(即插值函数)的图象并描绘出它与原始数据之间的位置关系如下图所示.第二节拟合(Fitting)正如上一节所讲插值主要用于获取一个给定的离散集合的局部逼近.如果要获取一个给定的离散集合的整体行为我们则采用拟合.也就是说给定一组离散数据xxxxnxyyyyny要求一个函数()xpy=按照最小二乘原理使得()iiniyxp=达到最小值这一方法被称为拟合法其中()xp被称为拟合曲线.值得注意的是拟合曲线()xpy=的确定可根据实验或实际情况来选择.通常我们考虑多项式拟合即选取拟合曲线()xp为m次多项式()()kkmkmxaxpxp==则问题变为求),,,,(mkakL=使得函数图光强度随时间的演化(),,===ikikmknimyxaaaafL达到最小值.确定),,,,(mkakL=的方法是让mjafj,,,,,L==于是得到mjxyaxnijiiknijkimk,,,,,L=====这是一个m元线性方程组.特别地当=m时即线性拟合此时拟合直线()axaxp=由系数a和a完全确定它们满足如下方程组=======niiniiniiiniiniiynaxayxxaxa通过引入向量矩阵=nxxxAMM和=nyyybM我们可以将该方程改写为如下的矩阵形式bAaaAATT=由此我们就可以得到拟合直线的系数为()bAAAaaTT=.下面是为求拟合直线的系数a和a所编写的MatLab函数文件(chapfitlinecm)functionyy=fitlinec(x,y)nx=length(x)ny=length(y)ifnx~=nywarning('Thelengthsofxandyshouldbesame')returnendn=min(nx,ny)ifn<error('ThenumberoftheDATAshouldbegreaterthan')returnendx=x(:n)y=y(:n)x=reshape(x,n,)y=reshape(y,n,)A=xones(n,)b=yB=A'*Ab=A'*byy=Bbyy=yy'和函数文件(chapfitlinepm)functionp=fitlinep(x,y,t)a=fitlinec(x,y)p=a*t一般地我们引入矩阵=nmnmnmnmmmxxxxxxxxxALMMOMMLL和=nyyybM则拟合曲线()xpm的系数,,,,aaaammL所满足的方程组可以改写为bAaaaAATmmT=M由此我们就可以得到拟合曲线的系数为()()bAAAaaaTTTmm=L.下面是为求拟合曲线()xpm的系数,,,,aaaammL所编写的MatLab函数文件(chapfitpolycm).functionyy=fitpolyc(x,y,m)nx=length(x)ny=length(y)ifnx~=nywarning('Thelengthsofxandyshouldbesame')returnendn=min(nx,ny)ifn<error('ThenumberoftheDATAshouldbegreaterthan')returnendifn<=mwarning('ThenumberoftheDATAshouldbegreaterthanm')returnendx=x(:n)y=y(:n)x=reshape(x,n,)y=reshape(y,n,)A=ones(n,)fori=:mA=x^iAendb=yB=A'*Ab=A'*byy=Bbyy=yy'和拟合曲线()xpm的函数文件(chapfitpolypm)functionpm=fitpolyp(x,y,t,m)a=fitpolyc(x,y,m)T=fori=:mT=t^iTendpm=a*T'特别地当=m时可以归结为直线拟合的情况.例某实验中测得一组数据其值如下x=y=.那么m次多项式拟合曲线与数据点的关系可以用如下程序的计算结果及图象来表述(chapfigurem)x=y=m=n=h=(max(x)min(x))nfori=:nt(i)=min(x)(i)*hpm(i)=fitpolyp(x,y,t(i),m)endplot(t,pm)holdonplot(x,y,'b*')运行结果给出拟合曲线与数据点的关系可以用如下图象来表述.图拟合曲线与数据点的关系第三节数值微分(NumericalDifferentiation)若函数()fx由表格给出求()fx在节点上的微商称这样的问题为数值微分即给定一组离散数据xxxxxnyyyyyn求函数()fx在节点()xini=,,,,L上的微商(其中()fxyii=).采用插值多项式来代替()fx去求函数在节点()xini=,,,,L处的微商.有三种思路:()先插值后计算微商()利用Taylor展开()插值与微商同时进行(如样条函数).下面我们利用Taylor展开来给出数值微分中的对称三点式、前向两点式、后向两点式和对称五点式的公式以及二阶导数的差分格式.已知xxhxhxxhxhyfffff求函数()fx在节点x处的微商(近似值).由()fx在x点处的Taylor展开()()()()()()()fxxfxxfxxfxxfx=′′′′′'!!LL知()()()()()()ffxhfhfxhfxhfxhfxh==′′′′′'()Ο()()()()()()()ffxhfhfxhfxhfxhfxh==′′′′′'()Ο()由此可以得到下面几个数值微分公式:()对称三点式由()得到()()()ffhfxhfxh=′′′'Ο即()()()′=′′′fxffhhfxhΟ取()′fxffh其截断误差为()Οh.()其中公式()被称为对称三点式公式.另外需要说明的是若()fx在xhxh,上是二次多项式则该估算式是精确的.实际上该估算式是xhxxh,,三点插值多项式的微商值.()前向两点式由()中f的表达式()()()()()ffhfxhfxhfxhfxh=′′′′′'()Ο知()()ffhfxh='Ο由此得到()′fxffh其截断误差为()Οh.()这里公式()被称为前向两点式公式.若()fx在xxh,上是一次多项式则该估算式是精确的.实际上该估算式是xxh,两点一次插值多项式的微商值.需要强调的是公式()的精确度比对称三点式低一个量级.()后向两点式由()中f的表达式()()()()()ffhfxhfxhfxhfxh=′′′′′'()Ο知()()ffhfxh='Ο由此得到()′fxffh其中截断误差为()Οh.()这里公式()被称为后向两点式公式.若()fx在xhx,上是一次多项式则该估算式是精确的.实际上该估算式是xhx,两点一次插值多项式的微商值.同样地公式()的精确度比对称三点式低一个量级与前向两点式是同一量级.()对称五点式为了得到更高的精确度同时利用()和()则有()()()ffhfxhfxh=′′′'Ο()()()ffhfxhfxh=′′′'Ο从中消去()′′′fx得到()()()()ffffhfxh='Ο由此得到()()′fxhffff()其中截断误差为()Οh我们称公式()被称为对称五点式公式.若()fx在xhxh,上是四次多项式则该估算式是精确的.实际上该估算式是xhxhxxhxh,,,,五点四次插值的微商值.注意该公式的精确度比对称三点式高两个量级.()二阶微商的数值微分利用()知()()()fffhfxhfxh=′′()Ο由此得到()′′fxfffh其中截断误差为()Οh.()()五点二阶微商的数值微分利用()和()知()()())(hxfhxfhfffΟ′′=()()())(hxfhxfhfffΟ′′=从中消去()xf)(我们得到()()()fffffhxf′′()其中截断误差为()hΟ.需要说明的是以上各公式均被称为差分公式.第四节数值积分(NumericalIntegration)本节目的是介绍计算积分()abfxdx的近似方法主要介绍梯形公式和Simpson公式.梯形公式利用定积分的几何意义可以导出如下梯形公式()()()()()()()abfxdxhfafahfahfanhfbOh=L()其中nabh=.抛物线(Simpson)公式计算积分()abfxdx的Simpson公式为()()()()()()hafhafhafhafafhdxxfba=()()()())()(hObfhnafhnafL其中nabh=.公式的推导:首先将区间ba,分成n等份每个小区间为,iixx,,,,=niL这样有()()dxxfdxxfiixxniba==.对于每个小区间,iixx上的积分()dxxfiixx我们用二次插值函数来代替其中的被积函数于是得到()()dxxfdxxfiixxniba==()()()()()=ixxiiixffxfxxii()()()()()hafhafhafhafafh=()()())()(bfhnafhnafL其中积分()dxxfiixx可以利用Maple来实现其中程序为:>restart:>l:=t>(t(xixi))*(txi)((xi(xixi))*(xixi)):>l:=t>(txi)*(txi)(((xixi)xi)*((xixi)xi)):>l:=t>(txi)*(t(xixi))((xixi)*(xi(xixi))):>simplify(int(l(t)*fil(t)*fil(t)*fi,t=xixi)):>factor()()fififi()xixi这里()()(),,===iixxiiixffffxffii.下面我们举几个利用MatLab计算数值积分的例子.例:计算积分dxex.显然该积分的精确值为=edxex.下面我们分别利用梯形公式和Simpson公式来计算.利用梯形公式的MatLab命令文件程序(命令文件chaptrapezoidm和函数文件integrandm):Usingtrapezoidruletocalculatetheintegrala=b=n=h=(ba)nint=integrand(a)fori=:nint=int*integrand(a(i)*h)endint=h*(intintegrand(b))formatlongint输出结果为int=其中被积函数的函数文件程序为:functiony=integrand(x)y=exp(x)利用Simpson公式的MatLab程序:(chapsimpsonm)UsingSimpsonruletocalcultetheintegrala=b=n=h=(ba)(*n)int=integrand(a)fori=:niint=int*integrand(a(*i)*h)endforj=:nint=int*integrand(a*j*h)endint=h*(intintegrand(b))int输出结果为int=.从上述的程序和计算结果可以看到利用梯形公式循环次所得的计算结果如果采用Simpson公式的话只需循环次即可.这说明Simpson公式具有比梯形公式更好地精度当然这也提高了计算的速度.例:计算积分()dxxx.这是一个瑕积分的例子有两个奇点=x和=x一般的处理方法是将积分区域分为若干个部分然后对每个积分作不同的变量替换消除掉奇点.就本例来说具体的处理形式为()dxxx()()()dxxxdxxxdxxx=()()()dxxxdxxxdxx=.作为练习.例:对于脉宽大于ps的光脉冲描述其在单模光纤内传输的方程是非线性薛定谔方程AATAAizAi||γβα=()其中A是脉冲包络的慢变振幅gvztT=是脉冲随gv移动的时间度量.方程()右端的三项分别对应于光脉冲在光纤中传输时的吸收效应、色散效应和非线性效应.根据入射脉冲的初始宽度T和峰值功率P决定脉冲在光纤内演变过程中是色散还是非线性效应起主要作用.我们首先将方程()无量纲化为此引入归一化振幅U和归一化初始脉宽T),(),(τταzUePzAz=TvztTTg==τ()式中P是入射脉冲的峰值功率指数因子代表光纤的损耗.这样方程()就变为UULeULzUiNLzD||)sgn(ατβ=()式中)sgn(=β它根据群速度色散参量β的符号来确定而||βTLD=PLNLγ=()分别被称为色散长度DL和非线性长度NLL它们分别表示沿光纤长L方向脉冲演化过程的长度量.根据色散长度DL、非线性长度NLL和光纤长度L的相对大小脉冲演变可分为下面讨论的四种不同的传输区域.()在DLL<<NLLL<<的情况当光纤长度DLL<<NLLL<<时色散和非线性效应都不起重要作用此时方程()右端的两项可被忽略.于是我们可以得到),(),(ττUzU=即脉冲在传输过程中保持其形状.在这个区域光纤不起太重要的作用只起传输光脉冲的作用除了由于吸收引起的脉冲能量的降低因而此区域对光通信系统是有益的.()在DLLNLLL<<的情况当光纤长度DLLNLLL<<时方程()右端的最后一项可被忽略此时在脉冲的演变过中群速度色散起主要作用而非线性效应相对较弱.由方程()可以看出当光纤和脉冲参量满足下述关系时适用于以色散为主的区域||<<=βγTPLLNLD()()在DLL<<NLLL的情况当光纤长度DLL<<NLLL时方程()右端的第一项较非线性项可以忽略.在这种情况下光纤中脉冲的演变过程自相位调制起主要作用它将导致脉冲的频谱展宽.由方程()可以看出当光纤和脉冲参量满足||>>=βγTPLLNLD()时适用于以非线性为主的区域.()在DLLNLLL的情况当光纤长度DLLNLLL时脉冲在光纤内传输过程中色散和非线性效应将共同起作用.群速度色散和自相位调制效应的互作用与群速度色散或自相位调制效应单独起作用相比较有不同的表现.在反常色散区(<β)光纤能维持光孤子.在正常色散区(>β)群速度色散和自相位调制可用来进行脉冲压缩.当群速度色散和自相位调制效应都重要时方程()对理解光纤中脉冲的演变是很有帮助的.在本例中我们将考虑在DLLNLLL<<的情况此时方程()或方程()右端的最后一项可被忽略可以写成TAzAi=β()方程()只包括了最低阶的群速度色散β尽管大多数实际情况下这一项的影响是主要的但有时也需要考虑更高阶的色散如三阶色散β.此时方程()可以改写为TAiTAzAi=ββ()此方程能利用傅立叶方法求解其解的形式为ωωωβωβωπdTiziziATzA=exp),(~),(()其中入射场),(TA的傅立叶变换),(~ωA为dTTiTAA)exp(),(),(~ωω=.()若入射场),(TA是确定的方程()就能用来研究高阶色散效应.下面我们来计算初始高斯脉冲的演化其中),(TA具有如下形式)exp(),(TTTA=.()在大多数情况下群速度色散β是主要的三阶色散β可以被忽略此时可以通过直接积分求得()的表达式.然而当=β时高阶色散效应β起主要作用此时我们很难求得积分()的表达式.这样为了看到高阶色散的作用数值计算积分()就显得非常重要.为此我们引入与高阶色散有关的三阶色散长度DL′||βTLD=′.()这样我们就可以将积分()和()写成如下形式=expexp),(~TTdTTTiTTTAωω()),(expexpTATdTiTωττωτ=和)(exp),(),(TdTTTiTzTiTATzAωωωβωπ=ϖϖϖϖπdTTiLziAD′=exp),(()其中()τϖττϖdiAexpexp),(=.()并假定>β.下面是我们编写的程序其中计算积分()的函数M文件为(chapgsintm)functionq=gsint(omega)b=n=ht=*bnt=b:ht:bint=forj=:nint=intht*(GS(t(j))*exp(i*omega*t(j))GS(t(j))*exp(i*omega*t(j)))endq=int而计算积分()的函数M文件为(chapgsintm)functionq=gsint(Z,t)b=n=hw=*bnw=b:hw:bint=forj=:nint=inthw*(gsint(w(j))*exp(i*w(j)^*Zi*w(j)*t)gsint(w(j))*exp(i*w(j)^*Zi*w(j)*t))endq=int*conj(int)(*pi)^最后可以得到运行程序(chapgsplotm)clearalla=n=hz=anb

职业精品

用户评论(1)

0/200
  • 泰瑞 2011-05-23 07:46:24

    不错

精彩专题

上传我的资料

热门资料

资料评价:

/35
仅支持在线阅读

意见
反馈

返回
顶部