《有限差分法基础课件.pptx》由会员分享，可在线阅读，更多相关《有限差分法基础课件.pptx（61页珍藏版）》请在三一办公上搜索。
1、第二章 有限差分法,主讲人：胡才博中国科学院大学地球科学学院中国科学院计算地球动力学重点实验室,第二章 有限差分法,2.1 有限差分法基础2.2 网格剖分2.3 差分格式2.4 差分方程2.5 应用实例,1. 地球内部介质，不仅存在纵向非均匀结构（一维地球模型），也存在横向非均匀结构（不同块体、断层系统）；2. 几何模型也呈现出相当的复杂性；3. 另外，边界条件和初始条件对于不同问题具有特殊性。,解析方法的局限性,Hu, C., Y. Cai, and Z. Wang (2012), Effects of large historical earthquakes, viscous relaxa
2、tion, and tectonic loading on the 2008 Wenchuan earthquake, Journal of Geophysical Research, 117, B06410, doi:10.1029/2011JB009046. (SCI, IF: 3.303),汶川大地震的动力学成因,对于存在复杂介质和几何、特殊边界条件和初始条件的实际地质问题，一般不存在解析解，需要近似的数值求解方法。有限差分方法是地球物理方法中最常见的一种。,有限差分方法（Finite Difference Method, FDM）是计算机数值模拟最早采用的方法，至今仍被广泛使用。,有限
3、差分方法的基本特点,该方法是一种直接将微分问题变为代数问题的近似数值解法，数学概念直观，表达简单，是发展较早且比较成熟的数值方法。,2.1 有限差分法基础,有限差分法以变量离散取值后对应的函数值来近似微分方程中独立变量的连续取值。 我们放弃了微分方程中独立变量可以取连续值的特征，而关注独立变量离散取值后对应的函数值。有限差分法的具体操作分为两个部分： (1)用差分代替微分方程中的微分，将连续变化的变量离散化，从而得到差分方程组的数学形式； (2)求解差分方程组。,有限差分方法的基本原理,该方法将求解域划分为差分网格，用有限个网格节点代替连续的求解域。有限差分方法以Taylor级数展开等方法，把
4、控制方程中的导数用网格节点上的函数值的差商代替进行离散，从而建立以网格节点上的函数值为未知数的代数方程组。,2.1 有限差分法基础,有限差分法的主要内容,建立地球物理问题的离散有限差分模型（1）如何根据问题的特点将定解区域做网格划分；（2）如何在所有网格节点上用有限差分格式对导数求近似，对函数、初始条件和边界条件求近似；（3）如何把原方程离散化为代数方程组，即有限差分方程组。,2.从理论上研究有限差分模型的形态，以保证计算过程的可行性和计算结果的正确性（1）解的相容性；（2）解的稳定性；（3）解的收敛性。,3. 如何数值求解差分方程组,网格剖分就是研究区域和边界的离散化1.矩形分割2.三角形分
5、割3.极网格分割,2.2 网格剖分,对地球物理问题的连续求解区域通过网格划分离散为空间上得一系列网格点，接下来需要利用一定的差分格式对偏微分方程组中的导数用差商进行近似，从而将偏微分方程组离散化为差分方程组。,对于函数f(x)，通常意义下的导数（微商）定义为：,2.3 差分格式,用Taylor级数展开可以给出微商的近似形式。,对于连续函数f(x)，它在相邻点上的值f(x+x)和f(x- x)可以用Taylor 级数展开为,dx变为x,如果x很小，f(x)可微，则以上级数收敛。次数越高，收敛级数的项的绝对值越小。,由（1）得到，,（1）,（2）,（3）,（4）,式中的O(x)项表示忽略掉的所有项
6、中的最大项的量级是x，也就是说，忽略掉这些项带来的误差中的最大项和x成正比。,由（4）给出导数的一阶精度（first order accurate）近似为：,（4）,（5）,（5）式称为向前差分格式（forward-difference formula）,由（2）式得到,由（7）式得到导数的另一个一阶精度近似：,（6）,（7）,（8）,（8）式称为向后差分形式（backward-difference formula）。,（1）式减去（2）式，得到：,（9）式中的O(x2)项表示忽略掉这些项带来的误差中的最大项和x2成正比。,由（9）式得到导数的二阶精度（second order accurat
7、e）近似为：,（10）式称为中心差分形式（central-difference formula）。,（9）,（10）,（1）式和（2）式相加，得到：,（12）式称为二阶导数的二阶精度中心差分形式。忽略x的四次方及更高阶项,(11),(12),f(xi+h)-f(xi): 节点xi的一阶向前差分f(xi)-f(xi-h): 节点xi的一阶向后差分f(xi+h)-f(xi-h): 节点xi的一阶中心差分前后是相对x轴正方向而言,总结：,1、向前差分形式：,2、向后差分形式：,3、中心差分形式：,单侧，一阶精度,单侧，一阶精度,对称，二阶精度,对于二阶导数,二阶精度,对一阶导数,定解问题的有限差分解
8、法1离散 x = ih, y= jh, i= 0, 1, 2,. n, h: 步长（正方形的边长）2根据泰勒级数建立差商格式：对于一维情况:在x处的一阶导数可以用,3. 建立和求解差商方程组,差分格式的另一种推导,为了寻求更精确的差分格式，我们引入两个待定常数，由泰勒展开，构造如下关系式,为了寻求更精确的差分格式，我们引入两个待定常数，由泰勒展开，构造如下关系式,回代（1）中，舍去高阶项,一阶偏导数中心差分的推导,（1）,（1）,二阶偏导数差分的推导,回代（1）中，舍去高阶项,二阶偏导数差分公式,一个例子：,等步长,一个例子：,局部节点离散化方程,总体节点离散化方程,总体节点离散化方程,f=0
9、时，变为泊松方程,f=q=0时，变为拉普拉斯方程,（1）第一类边界条件,（a） 直接转移法在图中网格是按正方形分割，步长为h。0点为靠近边界G的一个网格节点，1和2为边界节点。我们取最靠近0点的边界节点1上的函数值作为0点的函数值。即取01。这种方法称为直接转移法，又称为零次插值法。,边界条件的离散化的处理,狄利克莱问题,（b） 线性插值法先判断x方向的边界节点1和y方向的边界节点2哪一个更靠近0点。,如果1更靠近0点，则可以用x方向的线性插值给出0点的函数值,如果2更靠近0点，则可以用x方向的线性插值给出0点的函数值,（c） 双向插值法,变步长二次偏导数,（2） 第二类和第三类边界条件,对于
10、点O,过O点向边界G做垂线PQ交边界于Q，交网线段VR于P，OP=ah，PR=bh，VP=ch,因为P一般不是节点，其值应当以点和P、R点的插值给出,代入第二、三类边界条件,图中O与R重合,图中V与R点重合,（2） 第二类和第三类边界条件,2.4 差分方程,对于具体地球物理问题的偏微分方程组，利用上述差分格式，可以给出偏导数的微商近似，进一步得到差分方程组。,设f(P)是内部区域DI上定义的一个函数，设L(u)是一个微分算子，则以下表示了未知量u(P)的偏微分方程：,边界条件表示为以下方程：,设 和 分别表示区域D的北部节点和边界节点，则下式表示了以上偏微分方程的有限差分方程(finite-d
11、ifference equations, or finite-difference scheme)：,其中：以上差分方程是偏微分方程的有限差分近似，U是u的有限差分近似。差分方程要求U在所有节点上是u的很好的近似，并且方程所给出的有限差分近似解U是唯一的。,例子：对流方程（双曲型）的初值问题,差分方程,(13),(14),假定以上问题的解u(x,t)是充分光滑的，由Taylor级数展开有：,利用（15）和（17）式，得到：,如果u(x,t)是满足方程（13）的光滑解，则,代入（20），可以看出，偏微分方程（13）在(xj, tn)处可以近似地用下面的方程来代替：,其中是u(xj,tn)的近似值
12、。（21）式称为逼近方程（13）的有限差分方程，简称差分方程。用到的节点如图所示，,（21）式可以改写为便于计算的形式,(20),(21),(22),（22）式加上初始条件（14）的离散形式,(22),（23）,就可以按时间逐层推进，算出各层的值。这里的“层”表示在直线t=n上网格点的整体。差分方程（22）和初始条件的离散形式（23）结合在一起构成了一个差分格式。（22）给出了根据初始条件（23）来确定(j=0,1,)的一个算法，因此有时也称差分方程（22）为一个差分格式。差分格式包含了初始条件、边界条件的离散。,由第n个时间层推进到第n+1个时间层时，（22）提供了逐点直接计算n+1时的表达
13、式，因此称（22）为显式格式，在计算第n+1层时只用到了n层的数据，前后仅联系到连个时间层，因此（21）式为两层格式，更明确地称其为两层显式格式。,用（15）和（19），可以得到逼近方程（13）的另一差分格式,其中为网格比。此格式也是两层格式，称为中心差分格式。,用Taylor级数展开给出差分形式，用相应的差分形式逼近批未能微分方程（组），可以得到不同的差分格式，这一过程等价于用差商来近似微商得到相应的差分格式。,不同的差分格式的精度和误差不同。,例子：牛顿冷却定律：温度高于周围环境的物体向周围媒质传递热量逐渐冷却时所遵循的规律。当物体表面与周围存在温度差时，单位时间从单位面积散失的热量与温度
14、差成正比。,Tair,一阶常微分方程的数值解,如果不能简单地通过积分求解，则需要利用数值方法求解。首先对时间和温度进行离散：,利用向前差分形式：,得到以下的显式差分格式：,Tcap,利用该差分格式，我们可以计算任一时间和函数f的温度，但是随着计算的进行，误差O(dt)将会不断积累。,对于例子中的牛顿冷却问题，我们想知道液体的温度怎样随着时间和它与空气的温度差变化。,设温度差为,为冷却的时间尺度，这时,常微分方程为,初始条件：,以上的差分格式为：,该方程的解析解为：,该有限差分格式的近似程度如何？怎样选择t才能得到稳定解？,差分方程,该有限差分格式是否在t0时收敛？它是否和解析解得性质一样？,考
15、虑初始条件,得到,因此，对于第j个时刻,设t是到j时刻的总时间，,则：,对上式利用二项式展开,其中：,我们希望知道当dt0时差分格式的结果，这时相当于j,因此：,代入差分格式：,上式就是解析解,的Taylor展开。,结论：对于牛顿冷却问题，当时间步长t趋于零时，差分格式给出的数值解收敛于解析给出的严格解。,解析解是单调减小函数，数值解的性质怎样？保证数值解单调减小（Tj+1Tj）需要什么条件？,得到,保证Tj+1Tj的条件是：,因此，数值解只对一定取值范围的dt是单调减小的。,因此，数值结果是有条件稳定的。,差分方程,for i=1:nt+1 xt(i)=(i-1)*dt; T1(i)=t0*
16、exp(-xt(i)/tau);endplot(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-1Dclear;clc;figure(color,w);nt=8; % total time stepst0=1; % initial tempe
17、raturetau=0.7; % time constantdt=1.25; % time intervalT(1)=t0; for i=1:nt; xt(i)=(i-1)*dt;T(i+1)=T(i)-dt/tau*T(i);endxt(nt+1)=nt*dt;plot(xt,T,b.-);hold on,解析解,差分解,结果振荡但是收敛，前期结果不准确,结果收敛，前期结果不准确,结果收敛，误差逐渐缩小,结果收敛，结果基本吻合解析解,结果振荡，不收敛，计算结果错误,差分格式的性质分析,差分格式的性质分析,差分格式的性质分析,稳定性（stability）：如果偏微分方程的严格解析解有界，差分格
18、式给出的解也有界，称该差分格式是稳定的；如果差分格式给出的解是无界的，则称该差分格式是不稳定的。如果差分格式给出的解对于所有的时间步长和空间步长取值都是有界的，则称该差分格式是无条件稳定的；如果只是对时间步长和空间步长的部分取值有界，称它是有条件稳定的；如果对于所有的时间步长和空间步长取值都是无界的，则称差分格式是无条件非稳定的。稳定性反映了差分格式在计算中控制误差传递的能力，对偏微分方程有限差分方法研究具有重要意义。,例子：对流方程（双曲型）的初值问题,收敛性（convergence）：如果当时间和空间步长趋于零时，FDE解趋于PDE解，称该差分格式是收敛的。如果则称该差分格式是收敛的。对连
19、续形式的偏微分方程进行有限差分离散时，差分格式对最终计算结果有重要影响。收敛性表示当差分网格无限细化时，差分方程的解是否具有无限逼近偏微分方程的解的能力。收敛性的判断决定了一个差分格式能否被用来离散偏微分方程，不收敛的差分格式无疑是无法得到偏微分方程的真实解的。,利用向前差分，得到以下差分格式：,当时间步长dt趋近于零时，,差分格式的近似解趋近于方程的解析解,说明：该差分格式是收敛的。,收敛性实例,很重要的一点：相容性是差分格式的性质，它将差分格式和原始的偏微分方程联系在一起。稳定性和收敛性则是差分格式给出的数值解的性质。一般来讲，分析稳定性和相容性比较容易，证明收敛性有时是很困难的数学问题。
20、Lax等价定理（Lax equivalence theorem）：如果逼近一个给定问题的差分格式是相容的，那么该差分格式的收敛性与稳定性互为充分必要条件。该定理表明，对于一个具有相容性的差分格式，它的收敛性与稳定性是等价的。如果格式不稳定，则也不收敛。由于判断差分格式的收敛性相对比较困难，该定理提供了通过判断差分格式稳定性来确定收敛性的方法。因此，一般不再讨论收敛性问题，而是讨论差分格式的稳定性。该定理将收敛性与相容性和稳定性联系在一起，是很有用的一个定理。,差分格式的性质分析,拉普拉斯方程,U=0,U=0,U=100,U=0,拉普拉斯方程的有限差分法,问题：,2.5 应用实例,U=0,U=0
21、,U=100,U=0,2,1,4,3,6,5,8,7,10,9,12,11,14,13,15,组建A和B矩阵，求解线性方程组得到X,%Matlab 2Dclear;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;endfor i=1:7 a(15*i+15,15*i+14)=-0.25; a(15*i+15,15*i+30)=-0.25; a(15*i+
22、15,15*i)=-0.25; Enda(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; endfor i=122:134 a(i,i-1)=-0.25; a(i,i+1)=-0.25; a(i,i-15)=-0.25; endfor i=1:7for j=2:14;a(15*i
23、+j,15*i+j-1)=-0.25;a(15*i+j,15*i+j+1)=-0.25;a(15*i+j,15*i+j+15)=-0.25;a(15*i+j,15*i+j-15)=-0.25;endend,b=a(-1);c=zeros(135,1);for i=121:135 c(i,1)=25;endd=b*c;s=zeros(11,17);for i=2:16 s(11,i)=100; endfor i=1:9for j=1:15;s(i+1,j+1)=d(15*(i-1)+j,1);endend subplot(1,2,1),mesh(s)axis(0,17,0,11,0,100)su
24、bplot(1,2,2),contour(s,32),%Matlab 2Dclear;clc;figure(color,w); lx=17;ly=11; %v1=zeros(ly,lx); %for j=2:lx-1 v1(ly,j)=100;end %v2=v1;maxt=1;t=0;k=0;,while(maxt1e-6) %k=k+1maxt=0;for i=2:ly-1for j=2:lx-1;v2(i,j)=(v1(i,j+1)+v1(i+1,j)+v1(i-1,j)+v1(i,j-1)/4; %t=abs(v2(i,j)-v1(i,j);if(tmaxt) maxt=t;endendendv1=v2;end %,subplot(1,2,1),mesh(v1)axis(0,17,0,11,0,100)subplot(1,2,2),contour(v1,32),迭代解法:,K:迭代步数,K=419,2.5 应用实例,南加州一次未来大地震的强地面运动的数值模拟,盆地效应,Cui, 2013,Cui, 2013,Cui, 2013,Cui, 2013,总结：1、有限差分方法给出的数值解的精度取决于所用的差分形式（向前、向后、中心）。2、偏微分方程的显式有限差分格式通常是有条件稳定的，为了保证得到精确的数值解，最关键的是需要根据稳定性条件选取正确的空间和时间步长。,
