(整理)二维波动方程第一类吸收边界条件c++实现代码.

PML吸收边界条件影响因素分析丁科【摘要】影响完全匹配层方法吸收效果的主要因素有吸收层厚度和衰减系数,笔者通过地震波场数值模拟,讨论了吸收层厚度和衰减系数对吸收效果的影响.研究表明,在衰减系数一定时,吸收层厚度越大,边界反射越弱.吸收层的厚度一般在12～20道较为适当,在吸收层的厚度比较小时,随着衰减系数的增大,边界反射逐渐减弱甚至完全没有边界反射,但是若再进一步增大衰减系数则又会逐步出现边界反射.因此,在实际应用中应该注意衰减系数的选择,衰减系数一般在500～2 000之间,衰减系数较大反而影响其效果.%The main factors which influence the absorbing effects of the perfectly matched layer ( PML) method are the thickness of absorbing layer and the attenuation coefficient. Based on the numerical simulation of the seismic wave field, the author discusses the impact of the thickness of absorption layer and the attenuation coefficient on the absorption effect. The results show that, under the condition of constant attenuation coefficient, the larger the thickness of the absorption layer is, the weaker the boundary reflection becomes. The absorption layer thickness of 12 -20 traces is fairly appropriate. When the thickness of the absorption layer is smaller, the boundary reflection gets more and more weaker or even becomes zero with increasing attenuation coefficient. Nevertheless, if the attenuation coefficient is further enlarged, the boundary reflection will gradually appear again. Therefore, we must pay attention to the choice of attenuation coefficient in practice. The attenuation coefficient is generally between 500 and 2000. If theattenuation coefficient is larger than 3000, the absorption effect is not satisfactory.【期刊名称】《物探与化探》【年(卷),期】2012(036)004【总页数】5页(P623-627)【关键词】PML吸收边界;数值模拟;影响因素;波动方程【作者】丁科【作者单位】中南林业科技大学土木工程与力学学院,湖南长沙 410004【正文语种】中文【中图分类】P631.4波场数值模拟技术对于研究波在声介质、弹性介质、粘弹性介质等各种介质中的传播规律有着非常重要的作用，也是用来检验各种数据处理方法有效性的常用工具。

第七章 一维波动方程的傅氏解（20）一、内容摘要1．二阶线性偏微分方程可以分为如下四类：抛物型、双曲型、椭圆型和超双曲型方程。

抛物型：()2t xx yy zz u a u u u =++ 传导和扩散方程； 椭圆型：0xx yy zz u u u ++= Laplace 方程，稳态问题； 双曲型：()2tt xx yy zz u a u u u =++ 波动或弦振动方程。

2．一般地，要完全描写一个具有确定解的物理问题，在数学上就是要构成一个定解问题。

除了微分方程之外，构成定解问题还必须有边界条件和初始条件。

(1)初始条件：初始条件用于确定体系的历史状况，当所考察的物理现象 是随时间变化的时候，需要确定体系的初始条件来唯一确定地描述该现象。

(2)边界条件：体系的边界会影响体系的物理状态, 体系的边界情况由边界条件确定.边界条件反应体系和外界的界面上的情况。

常见的边界条件可以分为三类：①第一类边界条件：()(),,,|,r B u x y z t f r t ∈=r r． ②第二类边界条件：()|,r B u f r t n∈∂=∂r r． ③第三类边界条件：()()|,n r Bu cu f r t ∈+=r r． 上述三类边界条件,当函数(),0f r t =r时,分别称为第一、第二、第三类齐次边界条件。

3．定解问题问题的分类：数学物理方程（泛定方程）加上相应的定解条件一起构成了定解问题。

根据定解条件的不同，又可以把定解问题分为三类： 初值问题：定解条件仅有初值条件； 边值问题：定解条件仅有边值条件；混合问题：定界条件有初值条件也有边值条件。

4．分离变量法：（1）分离变量法的基本思想：将偏微分方程的问题转化为常微分方程的问题，先从中求出一些满足边界条件的特解，然后利用叠加原理，作出这些解的线性组合，令其满足余下的初始条件，从而得到定解问题的解。

（2）分离变量法的特点：把偏微分方程化为常微分方程,从而使问题的求解得以简化。

一类二维稳态对流——扩散方程的有限差分法一维稳态扩散方程描述了物质在一维空间中的扩散行为。

然而，在某些情况下，我们需要研究物质在二维平面中的扩散行为，例如热传导、流体传输等。

本文将介绍一类二维稳态对流-扩散方程的有限差分法。

二维稳态对流-扩散方程可以写作：∇·(D∇u) + ∇·(cu) + fu = 0 —— (1)其中，D是扩散系数，c是速度场，u是待求解的物理量，f是源项。

在这个方程中，第一项表示物质的扩散项，第二项表示对流项，第三项表示源项。

我们需要求解方程(1)，找到u的分布。

为了应用有限差分法来求解二维稳态对流-扩散方程，需要将二维空间离散化为一个网格。

假设我们将x方向离散为Nx个等距的节点，y方向离散为Ny个等距的节点，那么我们可以得到一个(Nx+1)×(Ny+1)的网格。

我们在网格节点上定义未知量u，然后将方程(1)对节点处的u进行离散化。

首先，我们对方程(1)的扩散项进行离散化。

我们使用五点差分格式来近似二维Laplace算符∇·(D∇u)。

对于网格节点(x,y)，我们可以得到以下差分格式：(Dij(xi+1,yj)ui+1,j + Dij(xi-1,yj)ui-1,j +Dij(xi,yj+1)ui,j+1 + Dij(xi,yj-1)ui,j-1 -4Dij(xi,yj)ui,j) / ∆x^2 + (Dij(xi,yj)ui,j) / ∆y^2其中，∆x和∆y是网格步长，Dij是扩散系数。

接下来，我们对方程(1)的对流项进行离散化。

我们使用中心差分格式来近似二维梯度算符∇·(cu)。

对于网格节点(x,y)，我们可以得到以下差分格式：(cxi+1/2,yj(ui+1,j - ui,j)) / ∆x + (cxi-1/2,yj(ui,j - ui-1,j)) / ∆x + (cyi,j+1/2(ui,j+1 - ui,j)) / ∆y + (cyi,j-1/2(ui,j - ui,j-1)) / ∆y其中，cxi+1/2,yj、cxi-1/2,yj、cyi,j+1/2和cyi,j-1/2是速度场在节点(x,y)处的中心点处的x和y分量。

数学物理方程小结 第七章 数学物理定解问题数学物理定解问题包含两个部分：数学物理方程（即泛定方程）和定解条件。

§7.1数学物理方程的导出一般方法： 第一确定所要研究的物理量u ,第二 分析体系中的任意一个小的部分与邻近部分的相互作用,根据物理规律, 抓住主要矛盾, 忽略次要矛盾。

（在数学上为忽略高级小量.）第三 然后再把物理量u 随时间，空间的变为通过数学算式表示出来, 此表示式即为数学物理方程。

（一） 三类典型的数学物理方程（1）波动方程： 0:),(:),(:22222222==∂∂-∂∂=∆-∂∂→f 当无外力时t x f xua t u 一维t r f u a tu 三维 此方程 适用于各类波动问题。

（特别是微小振动情况.）（2）输运方程： 0:).(:),(:2222==∂∂-∂∂=∆-∂∂→f 无外源时t x f xu a t u 一维t r f u a tu 三维 此方程 适用于热传导问题、扩散问题。

（3）Laplace 方程：.0(:0:).程时泊松方程退化拉氏方f f u 泊松方程u 拉氏方程t r ==∆=∆→稳定的温度和浓度分布适用的数学物理方程为Laplace 方程, 静电势u 在电荷密度为零处也满足Laplace 方程 。

§7.2定解条件定解条件包含初始条件与边界条件。

（1） 初始条件的个数等于方程中对时间最高次导数的次数。

例如波动方程应有二个初始条件, 一般选初始位移u （x,o ）和初始速度u t （x,0）。

而输运方程只有一个初始条件选为初始分布u （x,o ）,而Laplace 方程没有初始条件。

（2） 三类边界条件第一类边界条件： u( r ,t)|Σ = f （1） 第二类边界条件： u n |Σ = f （2） 第三类边界条件： ( u+Hu n )|Σ= f （3） 其中H 为常数. 7.3二阶线性偏微分方程分类判别式 ,,0,,0,,0221121222112122211212抛物型a a a 椭圆型a a a 双曲型a a a =-=∆<-=∆>-=∆ 波动方程是双曲型的,输运方程为抛物型的,而拉普拉斯方程为椭圆型的.7.4 达朗贝尔公式对一维无界的波动方程,当不考虑外力时,定解问题为()()()()()()()[]()⎰+-+++-====∂∂-∂∂atx at x t d aat x at x t x u 解为x x u x x u x u a t u ξξψϕϕψϕ2121,:0,0,022222对半无界问题作延拓处理:对第一类齐次边界条件作奇延拓,而对第二类齐次边界条件作偶延拓.第八章 分离变量法8.1分离变量法主要步骤:1.边界条件齐次化,对非齐次边界条件首先把它化为齐次的. •2.分离变量 u(x,t) =X(x) T(t) （1） [以后对三维问题也是如此]•3. 将（1）式代入原方程得出含任意常数λ的常微分方程, (称为本征方程) 而λ为本征值.•4.由齐次边界条件确定本征值,并求出本征方程.(得出的解为本征函数)•5.根据迭加原理把所有满足方程的线性无关解迭加后,就能得通解. •6.再由初始条件确定系数.一维波动方程在第一类齐次边界条件下的()()()()()()()()()4,sin 2:3,sin 22,sin 0,:1,sinsin cos ,:0011ξπξξψπξπξξϕϕππππd ln a n b 同样d ln l a x l xn a x u 代入边入边界l x n l at n b l at n a t x u 通解ln ln n n n n n ⎰⎰∑∑====⎪⎭⎫ ⎝⎛+=∞=∞=一维波动方程在第二类齐次边界条件下的通解:()()()()()()()()7.cos 2,cos 26.1,15,cossin cos .000000100ξπξξψπξπξξϕξξψξξϕπππd ln a n B d l n l A d l B d l A l x n l at n B l at n A t B A t x u ln ln ll n n n ⎰⎰⎰⎰∑====⎪⎭⎫ ⎝⎛+++=∞=一维输运方程在第一类齐次边界条件下的通解:()()()()9,sin 28,sin ,012⎰∑==⎪⎭⎫⎝⎛-∞=ln t l a n n n d ln l c lx n ec t x u ξπξξϕππ一维输运方程在第二类齐次边界条件下的通解:()()()()()11,cos 2,110,cos ,00002⎰⎰∑===⎪⎭⎫ ⎝⎛-∞=ln lt l a n n n d ln l c d l c lx n ec t x u ξπξξϕξξϕππ对其他的齐次边界条件,如本征函数已知也可直接求解,而对本征函数不熟则只能用分离变量法来求解. 8.2 非齐次边界条件的处理 常用方法有 1) 直线法 :对边界条件为: u(0,t)=g(t), u(L,t)=h(t) .令 ()()()()()x Lt g t h t g t x u t x v ---=,, ,可把边界条件化为齐次,但一般情况下方程变为非齐次. •只有当g,h 为常数时,方程才不变. 2) 特解法•把 u 化为两部分,令 u=v+w 使v 满足齐次边界条件与齐次方程,而使w 满足齐次方程与非齐次边界条件.下面通过实例来介绍此方法. • 例题 求解下列定解问题• U tt －a 2 U xx = 0 • U|x=0 =0, U|x=L = ASin ωt • U|t=0 = 0 , U t ∣t=0 = 0 •( 其中A 、ω为常数, 0＜x ＜L , 0＜ t ）•解：令 u=v+w ,使w 满足波动方程与非齐次边界条件,•得出()altaxA t x w ωωωsinsin sin,第九章 二阶常微分方程的级数解法 本征值问题一.拉普拉斯方程与亥姆霍斯方程在球坐标与柱坐标下分 离变量结果.1. 拉普拉斯方程在球坐标下的通解:()()()1,,1,,,1ϕϑϕϑim m l l L l l Y r B r A r u ∑⎪⎭⎫ ⎝⎛+=+其中Ylm为球函数,拉普拉斯方程在球坐标下的解不依赖于边界条件. 在轴对称时(1)式退化为()()()∑∞=+⎪⎭⎫⎝⎛+=012,cos ,l l l l l l P r B r A r u θθ2. 拉普拉斯方程在柱坐标下:()()()()()()()()()()()()()()()()()()..55.0:4,,0,ln :4;:3,04.01.3.022,1,0,sin cos 1.,,222222222''2程为m 阶Bessel方R m x dxdR x dx R d x 式为今x m F E R 式解为Bz A z Z 的解为R m d dR d R d Z Z m m m b m a z Z r R z u =-++==+=+===⎪⎪⎭⎫ ⎝⎛-++=-==+=ΦΦ=ρμρμρμρρρμλϕϕϕϕϕρ(5)式其解为m 阶Bessel 函数,解依赖于边界条件,当上下底为边界条件是齐次时, μ<0.对应的解是虚贝塞尔函数.3） 亥姆霍斯方程在球坐标与柱坐标下分离变量结果.在球坐标下:()()()ϕϑϕϑ,,,Y r R r u =其中Y 为球函数,R 为球贝塞尔函数.在柱坐标下:.()()()()()()()()()()()()()5.0:4,;4.01.3.022,1,0,sin cos 1.,,22222222222222''2=-++=-==⎪⎪⎭⎫ ⎝⎛--++=+==+=ΦΦ=R m x dxdR x dx R d x 式为今x k 令R m k d dR d R d Z Z m m m b m a z Z r R z u ρμνμρνρρρνλϕϕϕϕϕρ (5)式其解为m 阶Bessel 函数, 二、常微分方程的级数解法1. 掌握常点邻域的级数解法.2. 掌握正则奇点邻域的级数解法.3.知道无穷级数退化为多项式的方法. 三. 知道Sturm-Livouville 本征值问题的共同性质•当k(x),q(x)和ρ(x)都只取非负的值(≥0), Sturm-Livouville 方程共同性质为:•1)当k(x),k ’(x)和q(x)连续且x=a 和x=b 最多为一阶极点时,存在无限多个本征值及对应的本征函数:()()()()x y x y x y x y k k 321321,,≤≤≤≤≤λλλλ2)所有本征值λn ≥03)对应于不同本征值的本征函数带权正交()()()()m n dx x x y x y banm≠=⎰,0ρ4)本征函数族构成完备系()()∑∞==1n n n x y f x f第十章 球函数1. 轴对称的球函数当物理问题绕某一轴转动不变时,选此轴为z 轴这时物理量u 就与φ无关,m=0.此时球函数Y(θ,φ)就为L 阶勒让德多项式.即Y=P l (cos θ) 1) 勒让德多项式1. 勒让德多项式级数形式:()()()()()()1.!2!2!!22121202∑-=-----=l 或l n nl lnl x n l n l n n l x P 2. 勒让德多项式微分形式:()()()2.1!212l ll l l x dxd l x P -= 3.前几项为:P 0(x)= 1, P 1(x) =x=cos θ, •P 2(x)=(3x 2-1)/2, ….•一般勒让德多项式的幂次取决L•当L 为偶数时都为偶次幂项,L 为奇数时都为奇次幂项. 对特殊点x=1,0.()()()()()()()()(),!!2!!1210,00,1,11212n n P P x P x P P nn n l ll l --==-=-=-•4.勒让德多项式正交关系()lk l k l N dx x P x P δ211)(=⎰- (3) •5.勒让德多项式的模 122,1222+=+=l N l N l l (4) 6.广义傅里叶级数 :当f(x)在[-1,1]连续可导,且在x=-1与1有限时.()()()(),212111⎰∑-∞=+==dx x P x f l f x P f x f l l l l l (5) •7.在球坐标下Laplace 方程: △u= 0的通解为:轴对称()()()()()∑∑∑∞=+∞=-=+⎪⎭⎫⎝⎛+=⎪⎭⎫⎝⎛+=01017,cos 6,,l l l l l l l ll m lm l l l l P r B r A u Y r B r A r u θϕθθ (6)式有两系数需要两条件来确定,对球坐标有两自然边界条件,r=0与r →∞,球内解包含r=0,•u 有限, ()∑∞===0cos ,0l l ll l P r A u B θ (7)•而A l 由球面的边界条件确定,同样对球外区域两系数由球面的边界条件与r →∞, 两个条件确定. 8. 母函数()∑∞==+-02cos cos 211l l l P r r r θθ (8)9. 递推公式()()()()()()()0.12.2,112'1'1''1'111>-=+-+=++=+-+-++-l P P P l xP P P P x P l x lP x xP l l l l l l l l l l l二.连带勒让德函数•在一般情况下,物理量u 与φ有关,故球函数Y 是连带勒让德函数与周期函数的乘积. 1. 连带勒让德函数()[]()x P xm l m 221-=Θ (1)2.连带勒让德函数的微分表示()().1!21222lml m l lmml x dxd l x P --=++ (2) 从(2)可得当L 一定时,m 的取值为 m=0,1,2…L.共有L+1个值.而三角形式球函数Y （θ,φ）中,cosm φ,sinm φ为不同态,共有2L+1个态.3.正交关系()()()()()!!1223.2211m l m l l 模平方NN dx x P x P mllk ml m k m l -++==⎰-δ 4. 球函数Y 的两种表示形式. 第十一章 柱函数 一、 掌握三类柱函数的基本性质一般我们称Bessel 函数Jm(x)为第一类柱函数. 而把Neumann 函数Nm(x)称为第二类柱函数 . 1）对于第一类柱函数与第二类柱函数的线性组合.()()()()()()x iN x J x H x iN x J x H m m mm m m-=+=21称为第一种与第二种汉克尔函数.而汉克尔函数称为第三类柱函数2) x →0和x →∞时的行为()()()()()()()()()()()⎪⎭⎫⎝⎛---∞→⎪⎭⎫⎝⎛--∞→∞→∞→-→→→→==⎪⎭⎫ ⎝⎛--=⎪⎭⎫ ⎝⎛--=∞→∞→〉==4224210002lim ,2lim 42sin 2lim ,42cos 2lim lim ,lim 0.0lim ,1lim ππππππππππππm x i m x m x i m x m x m x m x m x m x x e xx H e xx H m x x x N m x x x J x J x N m x J x J3) 递推公式()()()()()()()[]()()()()()()()()()()()()4.3.212.1.211!21211!11'1'110122022x J xx J m x J x J x x J m x J 展开与把x J x x J x dxdxx J x k m k k x k m k dx d x J dx d m m m m m m m m m m mm k k k m k k kk m km m -+-+∞=-+∞=+=+-=-=-=⎪⎭⎫⎝⎛++Γ-=⎥⎥⎦⎤⎢⎢⎣⎡⎪⎭⎫ ⎝⎛++Γ-=⎥⎦⎤⎢⎣⎡∑∑4） 贝塞尔函数的零点对m 阶贝塞尔方程()()()()()()()()()()0)(::1.0.,0.00'222222====〉==-++ρμρμρμρμρμμρmm nm n m nmmJx 本征值x 记JJ R 件对柱侧面的齐次边界条时当x R m xdx dRxdx Rdx对第一类齐次边界条件 得出第n 个零点对第二类齐次边界条件 二．贝塞尔函数的正交关系 .• 对于不同本征值的同阶贝塞尔函数在区间 • [0,ρ0]上带权重ρ正交.• ()()()()()()1.][20nk m nm kmm nmNd J J δρρρμρμρ=⎰•• 2）广义傅里叶- 贝塞尔级数•()()()()()[]()()()()3.12.021ρρρμρρμρρd J f Nf J f f m nmm nn m n mn n ⎰∑==∞=• 3）Laplace 在柱坐标下的通解 • 轴对称m=0,柱内解为• 在侧面为第一类齐次边界条件时•()()()()()()()()()()2.,1.,101110000100⎪⎪⎭⎫ ⎝⎛⎥⎥⎦⎤⎢⎢⎣⎡⎪⎪⎭⎫ ⎝⎛+⎪⎪⎭⎫ ⎝⎛++=⎪⎪⎭⎫ ⎝⎛⎥⎥⎦⎤⎢⎢⎣⎡⎪⎪⎭⎫ ⎝⎛+⎪⎪⎭⎫ ⎝⎛=∑∑∞=∞=ρρρρR x J z R x sh B z Rx ch A z B A z u 条件时侧面为第二类齐次边界R x J z R x ch B z Rx sh A z u n n n n nn n n n n nn• 其中系数An,Bn 由上下底边界条件确定.• 在上下底为齐次边界条件时, μ≤ 0,R 的解为虚宗量贝塞尔函数.记为I m (x)• 同样可得Laplace 方程在柱内解 • 当轴对称时m=0• 上下底满足第一类齐次边界条件时解为•()()()()3.cos,:2.sin ,0001H z n H n I A z u 对第二类齐次边界条件H z n H n I A z u n n n n ππρρππρρ⎪⎭⎫⎝⎛=⎪⎭⎫⎝⎛=∑∑∞=∞=• 输运方程与波动方程在柱坐标下的解 • 1) 解的形式: u(r ,t)=T(t)v(r ) • V 满足亥姆霍兹方程.在侧面与上下底齐次边界条件下能完全确定本征值,例如上下底满足第一类齐次边界条件. 在轴对称情况下m=0 对输运方程柱内的解:上下底满足第一类齐次边界条件()()1.sin ,,2221,1000t H l x al n n nl n eH zl x J a t z u ⎥⎥⎦⎤⎢⎢⎣⎡⎪⎭⎫ ⎝⎛+⎪⎪⎭⎫ ⎝⎛-∞==∑⎪⎪⎭⎫ ⎝⎛=πρπρρρ波动方程在柱内的解:• 在上下底满足第一类齐次边界条件下•()[]()2002000000)(2.sin sin cos ,,⎪⎪⎭⎫ ⎝⎛+=⎪⎪⎭⎫ ⎝⎛+=∑∞ρπρρπρnnl n nlnl nl nl nl x H l k x J H z l at k b at k a t z u• 二维极坐标下的解: •• 侧面满足第一类齐次边界条件 ••()[]()∑∞=+=10000sin cos ,n n n n n n k J at k d at k c t u ρρ （3） •• 侧面满足第二类齐次边界条件•()[]()()4.sincos,1111ρρnnnnnnkJatkdatkct batu∑∞=+++=••第十二章积分变换法•一、傅里叶变换法••1。

第42卷第4期2003年12月石　油　物　探GEOPHYSICAL PROSPECTIN G FOR PETROL EUMVol.42,No.4Dec.,2003文章编号:100021441(2003)0420452204地震波方程人工边界的两种处理方法崔兴福1,2,张关泉2(1.中国石油辽河油田分公司勘探开发研究院计算所,辽宁盘锦124010;2.中国科学院数学与系统科学研究院计算数学与科学工程计算研究所,北京100080)摘要:在分析以往人工边界处理优缺点的基础上,提出了利用波动方程的频散关系式,得到可以吸收任何方向入射波的自适应校正吸收边界条件和旋转校正吸收边界条件。

同时,采用波阵面法和能流密度法判别入射波方向,克服了Pade近似吸收边界只对正入射波具有较好吸收性,而对非正入射的波吸收性不好的缺点。

数值试验结果表明了本方法的有效性。

关键词:自适应校正;旋转校正;波阵面;能流密度中图分类号:P63114 文献标识码:ATwo processing methods for artif icial boundary of seismic w ave equationCui Xingfu1,2,Zhang Guanquan2(puting Center of Exploration&Development Research Institute of CNPC Liaohe Oilfield Com pany,Panjin124010,China;2.Institute of Com putational Mathematics and Scientific/Engineering Computing,Academy of Math2 ematics and System Sciences,Chinese Academy of Sciences,Beijing100080,China)Abstract:In this paper,two absorbing boundary conditions,adaptive correction condition and rotation correction condi2 tion,were derived to absorb incident waves coming from any directions by using dispersion relation,based on an analy2 sis of the advantages and disadvantages of existing boundary conditions.The determination of incident wave direction by wave front and energy flux density was also described.Numerical ex periments were performed and their results were presented,which indicated the efficiency of these methods.K ey w ords:adaptive correction;rotation correction;wave front;energy flux density 实际人工模拟地震勘探是在半无界空间中进行的,而在计算机上进行数值模拟地震波在介质中的传播,只能在有限的模型空间中实现。