21有限差分法基础课堂

1、1,第二章,有限差分法,主讲人：胡才博,中国科学院大学地球科学学院,中国科学院计算地球动力学重点实,验室,2,第二章,有限差分法,2.1,有限差分法基础,2.2,网格剖分,2.3,差分格式,2.4,差分方程,2.5,应用实例,3,1,地球内部介质，不仅存在纵向非均匀结构（一维地球模型,也存在横向非均匀结构（不同块体、断层系统,2,几何模型也呈现出相当的复杂性,3,另外，边界条件和初始条件对于不同问题具有特殊性,解析方法的局限性,Hu, C,Y,Cai, and Z. Wang (2012), Effects of large historical earthquakes, viscous re

2、laxation, and tectonic loading,on the 2008 Wenchuan earthquake,Journal of Geophysical Research,117, B06410, doi:10.1029/2011JB009046,SCI, IF: 3.303,汶川大地震的动力学成因,4,对于存在复杂介质和几何、特殊边界条件和初始条件的实际地质问题,一般不存在解析解，需要近似的数值求解方法,有限差分方法是地球物理方法中最常见的一种,有限差分方法,Finite Difference Method, FDM,是计算机数值模拟最,早采用的方法，至今仍被广泛使用,有限

3、差分方法的基本特点,该方法是一种直接将微分问题变为代数问题的近似数值解法，数学,概念直观，表达简单，是发展较早且比较成熟的数值方法,2.1,有限差分法基础,5,有限差分法以变量离散取值后对应的函,数值来近似微分方程中独立变量的连续,取值,我们放弃了微分方程中独立变量可以取,连续值的特征，而关注独立变量离散取,有限差分方法的基本原理,该方法将求解域划分为差分网格，用有限个网格节点代替连续的求,解域。有限差分方法以,Taylor,级数展开等方法，把控制方程中的导数,用网格节点上的函数值的差商代替进行离散，从而建立以网格节点,上的函数值为未知数的代数方程组,2.1,有限差分法基础,6,有限差分法的主

4、要内容,1,建立地球物理问题的离散有限差分模型,1,如何根据问题的特点将定解区域做网格划分,2,如何在所有网格节点上用有限差分格式对导数求近似,对函数、初始条件和边界条件求近似,3,如何把原方程离散化为代数方程组，即有限差分方程组,2,从理论上研究有限差分模型的形态，以保证计算过程的可行性和计算结果的正确性,1,解的相容性,2,解的稳定性,3,解的收敛性,3,如何数值求解差分方程组,7,网格剖分就是研究区域和边界的离散化,1,矩形分割,2,三角形分割,3,极网格分割,2.2,网格剖分,8,对地球物理问题的连续求解区域通过网格划分离散为空间上得一系,列网格点，接下来需要利用一定的差分格式对偏微分

5、方程组中的导,数用差商进行近似，从而将偏微分方程组离散化为差分方程组,对于函数,f(x,通常意义下的导数（微商）定义为,lim,lim,lim,2,x,x,x,f,x,dx,f,x,dx,f,x,f,x,dx,dx,f,x,dx,f,x,dx,dx,当,dx0,时，以上三种形式都是微商的正确定义,如果,dx,是有限的，如何给出微商的近似定义,9,2.3,差分格式,用,Taylor,级数,展开可以给出微商的近似形式,对于连续函数,f(x,它在相邻点上的值,f(x,x,和,f(x,x,可以用,Taylor,级数展开为,dx,变为,x,2,2,3,3,4,4,2,3,4,2,2,3,3,4,4,2,

6、3,4,2,3,4,2,3,4,d,x,d,x,d,x,d,fx,x,fx,x,fx,fx,fx,fx,d,x,d,x,d,x,d,x,d,x,d,x,d,x,d,fx,x,fx,x,fx,fx,fx,fx,d,x,d,x,d,x,d,x,如果,x,很小,f(x,可微，则以上级数收敛,次数越高，收敛级数的项的绝对值越小,由,1,得到,2,2,3,3,4,4,2,3,4,2,3,4,d,x,d,x,d,x,d,fx,x,fx,x,fx,fx,fx,fx,d,x,d,x,d,x,d,x,d,fx,x,fx,fx,O,x,d,x,x,1,2,3,4,10,d,fx,x,fx,fx,O,x,d,x,x

7、,式中的,O(x,项表示忽略掉的所有项中的最大项的量级,是,x,也就是说，忽略掉这些项带来的误差中的最大,项和,x,成正比,d,f,x,x,f,x,f,x,d,x,x,由,4,给出导数的一阶精度,first order accurate,近似为,4,5,5,式称为,向前差分格式,forward-difference formula,由,2,式得到,2,2,3,3,4,4,2,3,4,2,3,4,d,x,d,x,d,x,d,fx,fxx,x,fx,fx,fx,fx,d,x,d,x,d,x,d,x,d,fx,fx,x,fx,O,x,d,x,x,由,7,式得到导数的另一个一阶精度近似,d,f,x,f

8、,x,x,f,x,d,x,x,6,7,8,8,式称为,向后差分形式,backward-difference formula,11,1,式减去,2,式，得到,2,2,d,fx,x,fxx,fx,O,x,d,x,x,9,式中的,O(x,2,项表示忽略掉这些项带来的误差中的最大项和,x,2,成正比,由,9,式得到导数的二阶精度,second order accurate,近似为,2,d,f,x,x,f,x,x,f,x,d,x,x,10,式称为,中心差分形式,central-difference formula,9,10,2,2,3,3,4,4,2,3,4,2,2,3,3,4,4,2,3,4,2,3,

9、4,2,3,4,d,x,d,x,d,x,d,fx,x,fx,x,fx,fx,fx,fx,d,x,d,x,d,x,d,x,d,x,d,x,d,x,d,fx,x,fx,x,fx,fx,fx,fx,d,x,d,x,d,x,d,x,1,2,12,2,2,3,3,4,4,2,3,4,2,2,3,3,4,4,2,3,4,2,3,4,2,3,4,d,x,d,x,d,x,d,fx,x,fx,x,fx,fx,fx,fx,d,x,d,x,d,x,d,x,d,x,d,x,d,x,d,fx,x,fx,x,fx,fx,fx,fx,d,x,d,x,d,x,d,x,1,2,2,2,2,2,2,d,fx,x,fx,fxx,f

10、x,O,x,d,x,x,1,式和,2,式相加，得到,2,2,2,2(,d,f,x,x,f,x,f,x,x,f,x,d,x,x,12,式称,为二阶导数的二阶精度中心差分形式,忽略,x,的四次方及更高阶项,11,12,f(x,i,h)-f(x,i,节点,x,i,的一阶,向前,差分,f(x,i,f(x,i,h,节点,x,i,的一阶,向后,差分,f(x,i,h)-f(x,i,h,节点,x,i,的一阶,中心,差分,前后是相对,x,轴,正方向,而言,13,d,f,x,x,f,x,f,x,d,x,x,d,f,x,f,x,x,f,x,d,x,x,2,2,2,2(,d,f,x,x,f,x,f,x,x,f,x,d

11、,x,x,2,d,f,x,x,f,x,x,f,x,d,x,x,总结,1,向前差分形式,2,向后差分形式,3,中心差分形式,单侧，一阶精度,单侧，一阶精度,对称，二阶精度,对于二阶导数,二阶精度,对一阶导数,14,定解问题的有限差分解法,1,离散,x = ih, y= jh, i= 0,1,2,n, h,步长（正方形的边长,2,根据泰勒级数建立差商格式：对于一,维情况,在,x,处的一阶导数可以用,3,建立和求解差商方程组,15,差分格式的另一种推导,为了寻求更精确的差分格式，我们引入两个,待定常数，由泰勒展开，构造如下关系式,i,j,i+1,j,i-1,j,0,1,2,3,4,i,j,i+1,j

12、-1,i-1,j-1,i,j+1,i+1,j+1,i-1,j+1,i-1,i,i+1,j-1,j,j+1,h,1,h,3,h,2,h,4,16,为了寻求更精确的差分格式，我们引入两个,待定常数，由泰勒展开，构造如下关系式,回代,1,中，舍去高阶项,一阶偏导数中心差分的推导,1,17,1,二阶偏导数差分的推导,回代,1,中，舍去高阶项,二阶偏导数差分公式,18,一个例子,等步长,19,一个例子,局部节点离散化方程,总体节点离散化方程,总体节点离散化方程,f=0,时，变为泊松方程,f=q=0,时，变为拉普拉斯方程,20,1,第一类边界条件,a,直接转移法,在图中网格是按正方形分割,步长为,h,0,

13、点为靠近边界,G,的一个网格节点,1,和,2,为边界节点。我们取最靠近,0,点的,边界节点,1,上的函数值作为,0,点的函数值,即取,0,1,这种方法称为直接转移法，又称为零次插值法,边界条件的离散化的处理,G,g,s,狄利克莱问题,21,b,线性插值法,先判断,x,方向的边界节点,1,和,y,方向的边界节点,2,哪一个更靠近,0,点,如果,1,更靠近,0,点，则可以用,x,方向的线性插,值给出,0,点的函数值,如果,2,更靠近,0,点，则可以用,x,方向的线性插,值给出,0,点的函数值,22,c,双向插值法,i,j,i+1,j,i-1,j,0,1,2,3,4,i,j,i+1,j-1,i-1,

14、j-1,i,j+1,i+1,j+1,i-1,j+1,i-1,i,i+1,j-1,j,j+1,h,1,h,3,h,2,h,4,变步长二次偏导数,23,2,第二类和第三类边界条件,对于点,O,过,O,点向边界,G,做垂线,PQ,交边界于,Q,交网线段,VR,于,P,OP=ah,PR=bh,VP=ch,因为,P,一般不是节点，其值应当以点和,P,R,点的,插值给出,代入第二、三类边界条件,24,图中,O,与,R,重合,图中,V,与,R,点重合,2,第二类和第三类边界条件,25,2.4,差分方程,对于具体地球物理问题的偏微分方程组，利用上述差分格式，可以给出偏,导数的微商近似，进一步得到差分方程组,设

15、,f(P,是内部区域,D,I,上定义的一个函数，设,L(u,是一个微分算子，则以下表,示了未知量,u(P,的偏微分方程,边界条件表示为以下方程,设,和,分别表示区域,D,的北部节点和边界节点，则下式表示了以,上偏微分方程的有限差分方程,finite-difference equations, or finite,difference scheme,I,L,u,P,f,P,PD,B,B,u,P,g,P,PD,I,LU,f,P,P,D,B,B,U,g,P,PD,其中：以上差分方程是偏微分方程的有限差分近似,U,是,u,的有限差分近似,差分方程要求,U,在所有节点上是,u,的很好的近似，并且方程所给

16、出的有限差,分近似解,U,是唯一的,I,D,B,D,26,例子：对流方程（双曲型）的初值问题,差分方程,0,0,0,u,u,a,x,R,t,t,x,u,x,g,x,x,R,13,14,假定以上问题的解,u(x,t,是充分光滑的，由,Taylor,级数展开有,27,利用,15,和,17,式，得到,如果,u(x,t,是满足方程,13,的光滑解，则,代入,20,，可以看出，偏微分方程,13,在,x,j,t,n,处可以近似地,用下面的方程来代替,其中是,u(x,j,t,n,的近似值。,21,式称为逼近方程,13,的有限差分方程，简,称差分方程。用到的节点如图所示,21,式可以改写为便于计算的形式,1,

17、1,n,j,n,j,n,j,n,j,n,j,u,x,t,u,x,t,u,xt,u,x,t,u,u,a,a,O,h,h,t,x,20,0,n,j,u,u,a,t,x,1,1,0,0,1,2,n,0,1,2,n,n,n,n,j,j,j,j,u,u,u,u,a,h,j,21,1,1,n,n,n,n,j,j,j,j,u,u,a,u,u,22,28,22,式加上初始条件,14,的离散形式,0,0,1,2,j,j,u,j,1,1,n,n,n,n,j,j,j,j,u,u,a,u,u,22,23,就可以按时间逐层推进，算出各层的值。这里的“层”表示在直线,t=n,上网格,点的整体。差分方程,22,和初始条件的

18、离散形式,23,结合在一起构成了一,个差分格式。,22,给出了根据初始条件,23,来确定,j=0,1,的一个算法,因此有时也称差分方程,22,为一个差分格式。差分格式包含了初始条件、边,界条件的离散,由第,n,个时间层推进到第,n+1,个时间层时，,22,提供了逐点直接计算,n+1,时的,表达式，因此称,22,为显式格式，在计算第,n+1,层时只用到了,n,层的数据，前,后仅联系到连个时间层，因此,21,式为两层格式，更明确地称其为两层显式,格式,用,15,和,19,，可以得到逼近方程,13,的另一差分格式,1,1,1,0,2,0,1,2,n,0,1,2,n,n,n,n,j,j,j,j,u,u

19、,u,u,a,h,j,1,1,1,2,n,n,n,n,j,j,j,j,a,u,u,u,u,其中,为网格比。此格式也是两层格式，称为中心差分格式,29,用,Taylor,级数展开给出差分形式，用相应的差分形式逼近批未能微分方,程（组），可以得到不同的差分格式，这一过程等价于用差商来近似微,商得到相应的差分格式,不同的差分格式的精度和误差不同,30,例子：牛顿冷却定律,温度高于周围环境的物体向周围媒质传递热量,逐渐冷却时所遵循的规律。当物体表面与周围存在温度差时，单位时,间从单位面积散失的热量与温度差成正比,T,air,一阶常微分方程的数值解,dT,f,T,t,dt,如果不能简单地通过积分求解，则

20、需要利用数值方法求解。首先对时间和温度进,行离散,0,j,j,j,t,t,j,t,T,T,t,利用向前差分形式,1,j,j,j,t,t,T,T,d,T,O,t,d,t,t,得到以下的显式差分格式,1,j,j,j,j,T,T,t,f,Tt,T,cap,31,利用该差分格式，我们可以计算任一时间和函数,f,的温度,但是随着计算的进行，误差,O(dt,将会不断积累,1,j,j,j,j,T,T,t,f,Tt,对于例子中的牛顿冷却问题，我们想知道液体的温度怎,样随着时间和它与空气的温度差变化,设温度差为,c,a,p,a,i,r,T,T,T,为冷却的时间尺度，这时,f,Tt,T,常微分方程为,d,T,T,

21、d,t,初始条件,0,0,T,T,t,以上的差分格式为,1,1,j,j,j,j,t,t,T,T,T,T,该方程的解析解为,0,e,x,p,T,t,T,t,该有限差分格式的近似程度如何？怎样选择,t,才能得到稳定解,32,差分方程,1,1,j,j,t,T,T,该有限差分格式是否在,t0,时收敛？它是否和解析解得性质一样,考虑初始条件,0,0,T,T,t,得到,1,0,1,t,T,T,2,2,1,0,1,1,t,t,T,T,T,因此，对于第,j,个时刻,0,1,j,j,t,T,T,设,t,是到,j,时刻的总时间,t,t,j,则,0,1,j,j,t,T,T,j,1,2,2,0,11,1,1,2,j,

22、j,j,j,j,j,t,t,TT,j,j,对上式利用二项式展开,其中,j,j,r,j,r,r,我们希望知道当,dt0,时差分格式的结果，这时相当于,j,1,2,1,r,j,jj,j,j,r,j,j,r,33,因此,r,j,j,r,r,代入差分格式,2,1,2,0,2,0,1,1,2,1,1,2,j,j,j,t,j,t,T,T,j,j,t,t,T,T,上式就是解析解,0,e,x,p,T,t,T,t,的,Taylor,展开,结论：对于牛顿冷却问题，当时间步长,t,趋于零时，差分格式给出的,数值解收敛于解析给出的严格解,解析解是单调减小函数，数值解的性质怎样？保证数值解单调减小,T,j+1,T,j,

23、需要什么条件,1,2,2,0,11,1,1,2,j,j,j,j,j,j,t,t,TT,j,j,得到,j,j,r,j,r,r,1,2,1,r,j,jj,j,j,r,j,j,r,34,1,1,j,j,t,T,T,0,1,1,t,0,t,保证,T,j+1,T,j,的条件是,因此，数值解只对一定取值范围的,dt,是单调减小的,如果,2,t,1,1,0,t,1,1,t,数值解振荡，但是结果会收敛,如果,2,t,1,1,t,1,1,t,数值解振荡并且发散,因此，数值结果是有条件稳定的,35,差分方程,for i=1:nt+1,xt(i)=(i-1)*dt,T1(i)=t0*exp(-xt(i)/tau,e

24、nd,plot(xt,T1,r.-);hold on,set(gca,DataAspectRatio,(max(xt),min(xt)/(max(T)-min(T)/3 1 1,xlabel(Time (s),Fontname,times new,roman,FontSize,14,ylabel(Temperature,Fontname,times new,roman,FontSize,14,Malab-1D,clear,clc,figure(color,w,nt=8; % total time steps,t0=1; % initial temperature,tau=0.7; % time

25、 constant,dt=1.25; % time interval,T(1)=t0,for i=1:nt,xt(i)=(i-1)*dt,T(i+1)=T(i)-dt/tau*T(i,end,xt(nt+1)=nt*dt,plot(xt,T,b.-);hold on,d,T,T,d,t,0,e,x,p,T,t,T,t,1,1,j,j,t,T,T,解析解,差分解,36,结果振荡但是收敛,前期结果不准确,结果收敛,前期结果不准确,37,结果收敛,误差逐渐缩小,结果收敛,结果基本吻合解析解,38,结果振荡，不收敛，计算结果错误,39,差分格式的性质分析,40,41,42,差分格式的性质分析,43,差

26、分格式的性质分析,稳定性,stability,：如果偏微分方程的严格解析解有界，差分格式给出,的解也有界，称该差分格式是稳定的；如果差分格式给出的解是无界的,则称该差分格式是不稳定的,如果差分格式给出的解对于所有的时间步长和空间步长取值都是有界的,则称该差分格式是无条件稳定的；如果只是对时间步长和空间步长的部,分取值有界，称它是有条件稳定的；如果对于所有的时间步长和空间步,长取值都是无界的，则称差分格式是无条件非稳定的,稳定性反映了差分格式在计算中控制误差传递的能力，对偏微分方程有,限差分方法研究具有,重要意义,0,0,0,u,u,a,x,t,t,x,u,x,g,x,x,R,R,例子：对流方程

27、（双曲型）的初值问题,44,0,0,0,u,u,a,x,t,t,x,u,x,g,x,x,R,R,45,收敛性,convergence,如果当时间和空间步长趋于零时,FDE,解趋于,PDE,解，称该差分格式是收敛的,如果,则称该差分格式是收敛的,对连续形式的偏微分方程进行有限差分离散时，差分格式,对最终计算结果有重要影响。收敛性表示当差分网格无限,细化时，差分方程的解是否具有无限逼近偏微分方程的解,的能力。收敛性的判断决定了一个差分格式能否被用来离,散偏微分方程，不收敛的差分格式无疑是无法得到偏微分,方程的真实解的,0,0,l,i,m,0,m,m,i,jk,i,jk,h,t,U,u,46,d,T

28、,T,d,t,2,1,2,0,2,0,1,1,2,1,1,2,j,j,j,t,j,t,T,T,j,j,t,t,T,T,利用向前差分，得到以下差分格式,1,1,j,j,j,j,t,t,T,T,T,T,当时间步长,dt,趋近于零时,差分格式的近似解趋近于方程的解析解,0,e,x,p,T,t,T,t,说明：该差分格式是收敛的,收敛性实例,47,很重要的一点：相容性是差分格式的性质，它将差分格式和原始的偏微分方程,联系在一起。稳定性和收敛性则是差分格式给出的数值解的性质,一般来讲，分析稳定性和相容性比较容易，证明收敛性有时是很困难的数学问,题,Lax,等价定理,Lax equivalence theo

29、rem,：如果逼近一个给定问题的差分格式是,相容的，那么该差分格式的收敛性与稳定性互为充分必要条件,该定理表明，对于一个具有相容性的差分格式，它的收敛性与稳定性是等价的,如果格式不稳定，则也不收敛。由于判断差分格式的收敛性相对比较困难，该,定理提供了通过判断差分格式稳定性来确定收敛性的方法。因此，一般不再讨,论收敛性问题，而是讨论差分格式的稳定性,该定理将收敛性与相容性和稳定性联系在一起，是很有用的一个定理,差分格式的性质分析,48,1,1,1,1,4,U,ij,U,ij,U,i,j,U,i,j,U,i,j,2,0,U,拉普拉斯方程,0,2,5,1,0,2,5,1,0,2,5,1,0,2,5,

30、1,0,U,i,j,U,i,j,U,i,j,U,i,j,U,i,j,U=0,U=0,U=100,U=0,2,0,U,拉普拉斯方程的,有限差分法,问题,2.5,应用实例,49,U=0,U=0,U=100,U=0,2,0,U,2,1,4,3,6,5,8,7,10,9,12,11,14,13,15,0,2,5,1,0,2,5,1,0,2,5,1,0,2,5,1,0,U,i,j,U,i,j,U,i,j,U,i,j,U,i,j,50,A,X,B,1,3,5,1,3,5,A,A,1,3,5,1,X,X,1,3,5,1,B,B,组建,A,和,B,矩阵，求解线性方程组得到,X,51,Matlab 2D,cle

31、ar,clc,figure(color,w,a=zeros(135,135,for i=1:135,a(i,i)=1,end,for i=1:7,a(15*i+1,15*i+2)=-0.25,a(15*i+1,15*i+16)=-0.25,a(15*i+1,15*i-14)=-0.25,end,for i=1:7,a(15*i+15,15*i+14)=-0.25,a(15*i+15,15*i+30)=-0.25,a(15*i+15,15*i)=-0.25,End,a(1,2)=-0.25,a(1,16)=-0.25,a(121,122)=-0.25,a(121,106)=-0.25,a(135,134)=-0.25,a(135,120)=-0.25,a(15,14)=-0.25,a(15,30)=-0.25,for i=2:14,a(i,i-1)=-0.25,a(i,i+1)=-0.25,a(i,i+15)=-0.25,end,for i=122:134,a(i,i-1)=-0.25,a(i,i+1)=-0.25,a(i,i-15)=-0.25,end,for i=1:7,for j=2:14,a(