2.1有限差分法基础4

2021/6/161第二章 有限差分法主讲人：胡才博主讲人：胡才博中国科学院大学地球科学学院中国科学院大学地球科学学院中国科学院计算地球动力学重点实中国科学院计算地球动力学重点实验室验室2021/6/162第二章 有限差分法 2.1 有限差分法基础 2.2 网格剖分 2.3 差分格式 2.4 差分方程 2.5 应用实例2021/6/1631.地球内部介质，不仅存在纵向非均匀结构（一维地球模型），也存在横向非均匀结构（不同块体、断层系统）；2.几何模型也呈现出相当的复杂性；3.另外，边界条件和初始条件对于不同问题具有特殊性。解析方法的局限性解析方法的局限性Hu,C.,Y.Cai,and Z.Wang(2012),Effects of large historical earthquakes,viscous relaxation,and tectonic loading on the 2008 Wenchuan earthquake,Journal of Geophysical Research,117,B06410,doi:10.1029/2011JB009046.(SCI,IF:3.303)汶川大地震的动力学成因2021/6/164对于存在复杂介质和几何、特殊边界条件和初始条件的实际地质问题，一般不存在解析解，需要近似的数值求解方法。有限差分方法是地球物理方法中最常见的一种。有限差分方法（Finite Difference Method,FDM）是计算机数值模拟最早采用的方法，至今仍被广泛使用。有限差分方法的基本特点该方法是一种直接将微分问题变为代数问题的近似数值解法，数学概念直观，表达简单，是发展较早且比较成熟的数值方法。2.1 有限差分法基础2021/6/165 有限差分法以变量离散取值后对应的函数值来近似微分方程中独立变量的连续取值。我们放弃了微分方程中独立变量可以取连续值的特征，而关注独立变量离散取值后对应的函数值。有限差分法的具体操作分为两个部分：(1)用差分代替微分方程中的微分，将连续变化的变量离散化，从而得到差分方程组的数学形式；(2)求解差分方程组。有限差分方法的基本原理该方法将求解域划分为差分网格，用有限个网格节点代替连续的求解域。有限差分方法以Taylor级数展开等方法，把控制方程中的导数用网格节点上的函数值的差商代替进行离散，从而建立以网格节点上的函数值为未知数的代数方程组。2.1 有限差分法基础2021/6/166有限差分法的主要内容建立地球物理问题的离散有限差分模型（1）如何根据问题的特点将定解区域做网格划分；（2）如何在所有网格节点上用有限差分格式对导数求近似，对函数、初始条件和边界条件求近似；（3）如何把原方程离散化为代数方程组，即有限差分方程组。2.从理论上研究有限差分模型的形态，以保证计算过程的可行性和计算结果的正确性（1）解的相容性；（2）解的稳定性；（3）解的收敛性。3.如何数值求解差分方程组2021/6/167 网格剖分就是研究区域和边界的离散化 1.矩形分割 2.三角形分割 3.极网格分割2.2 网格剖分2021/6/168对地球物理问题的连续求解区域通过网格划分离散为空间上得一系列网格点，接下来需要利用一定的差分格式对偏微分方程组中的导数用差商进行近似，从而将偏微分方程组离散化为差分方程组。对于函数f(x)，通常意义下的导数（微商）定义为：()()lim()()lim()()lim2xxxf xdxf xdxf xf xdxdxf xdxf xdxdx 当dx0时，以上三种形式都是微商的正确定义。如果dx是有限的，如何给出微商的近似定义？2021/6/1692.3 差分格式用Taylor级数级数展开可以给出微商的近似形式。对于连续函数f(x)，它在相邻点上的值f(x+x)和f(x-x)可以用Taylor 级数展开为dx变为x223344234223344234()()()()()()()()()2!3!4!()()()()()()()()()2!3!4!dxdxdxdf xxf xxf xf xf xf xdxdxdxdxdxdxdxdf xxf xxf xf xf xf xdxdxdxdx如果x很小，f(x)可微，则以上级数收敛。次数越高，收敛级数的项的绝对值越小。由（1）得到，223344234()()()()()()()()()2!3!4!()()()()dxdxdxdf xxf xxf xf xf xf xdxdxdxdxdf xxf xf xOxdxx（1）（2）（3）（4）2021/6/1610()()()()df xxf xf xOxdxx式中的O(x)项表示忽略掉的所有项中的最大项的量级是x，也就是说，忽略掉这些项带来的误差中的最大项和x成正比。()()()df xxf xf xdxx由（4）给出导数的一阶精度（first order accurate）近似为：（4）（5）（5）式称为向前差分格式向前差分格式（forward-difference formula）由（2）式得到223344234()()()()(-)()-()()-()2!3!4!()()()+()dxdxdxdf xf xxxf xf xf xf xdxdxdxdxdf xf xxf xOxdxx 由（7）式得到导数的另一个一阶精度近似：()(-)()df xf xxf xdxx（6）（7）（8）（8）式称为向后差分形式向后差分形式（backward-difference formula）。2021/6/1611（1）式减去（2）式，得到：2()(-)()=+()2df xxf xxf xOxdxx（9）式中的O(x2)项表示忽略掉这些项带来的误差中的最大项和x2成正比。由（9）式得到导数的二阶精度（second order accurate）近似为：()(-)()2df xxf xxf xdxx（10）式称为中心差分形式中心差分形式（central-difference formula）。（9）（10）223344234223344234()()()()()()()()()2!3!4!()()()()()()()()()2!3!4!dxdxdxdf xxf xxf xf xf xf xdxdxdxdxdxdxdxdf xxf xxf xf xf xf xdxdxdxdx（1）（2）2021/6/1612223344234223344234()()()()()()()()()2!3!4!()()()()()()()()()2!3!4!dxdxdxdf xxf xxf xf xf xf xdxdxdxdxdxdxdxdf xxf xxf xf xf xf xdxdxdxdx（1）（2）2222()2()(-)()=+()()df xxf xf xxf xOxdxx（1）式和（2）式相加，得到：222()2()(-)()()df xxf xf xxf xdxx（12）式称为二阶导数的二阶精度中心差分形式为二阶导数的二阶精度中心差分形式。忽略x的四次方及更高阶项(11)(12)f(xi+h)-f(xi):节点xi的一阶向前差分 f(xi)-f(xi-h):节点xi的一阶向后差分 f(xi+h)-f(xi-h):节点xi的一阶中心差分 前后是相对x轴正方向而言2021/6/1613()()()df xxf xf xdxx()(-)()df xf xxf xdxx222()2()(-)()()df xxf xf xxf xdxx()(-)()2df xxf xxf xdxx总结：总结：1、向前差分形式：2、向后差分形式：3、中心差分形式：单侧，一阶精度单侧，一阶精度对称，二阶精度对于二阶导数二阶精度对一阶导数2021/6/1614 定解问题的有限差分解法1离散 x=ih,y=jh,i=0,1,2,.n,h:步长（正方形的边长）2根据泰勒级数建立差商格式：对于一维情况:在x处的一阶导数可以用3.建立和求解差商方程组2021/6/1615差分格式差分格式的另一种推导的另一种推导为了寻求更精确的差分格式，我们引入两个待定常数，由泰勒展开，构造如下关系式(i,j)(i+1,j)(i-1,j)01234(i,j)(i+1,j-1)(i-1,j-1)(i,j+1)(i+1,j+1)(i-1,j+1)i-1ii+1j-1jj+1h1h3h2h42021/6/1616为了寻求更精确的差分格式，我们引入两个待定常数，由泰勒展开，构造如下关系式回代（1）中，舍去高阶项一阶偏导数中心差分的推导一阶偏导数中心差分的推导（1）2021/6/1617（1）二阶偏导数差分的推导二阶偏导数差分的推导回代（1）中，舍去高阶项二阶偏导数差分公式2021/6/1618一个例子：等步长2021/6/1619一个例子：局部节点离散化方程总体节点离散化方程总体节点离散化方程f=0时，变为泊松方程f=q=0时，变为拉普拉斯方程2021/6/1620（1）第一类边界条件（a）直接转移法在图中网格是按正方形分割，步长为h。0点为靠近边界G的一个网格节点，1和2为边界节点。我们取最靠近0点的边界节点1上的函数值作为0点的函数值。即取01。这种方法称为直接转移法，又称为零次插值法。边界条件的离散化的处理边界条件的离散化的处理()Gg s狄利克莱问题2021/6/1621（b）线性插值法先判断x方向的边界节点1和y方向的边界节点2哪一个更靠近0点。如果1更靠近0点，则可以用x方向的线性插值给出0点的函数值如果2更靠近0点，则可以用x方向的线性插值给出0点的函数值2021/6/1622（c）双向插值法(i,j)(i+1,j)(i-1,j)01234(i,j)(i+1,j-1)(i-1,j-1)(i,j+1)(i+1,j+1)(i-1,j+1)i-1ii+1j-1jj+1h1h3h2h4变步长二次偏导数2021/6/1623（2）第二类和第三类边界条件 对于点O过O点向边界G做垂线PQ交边界于Q，交网线段VR于P，OP=ah，PR=bh，VP=ch因为P一般不是节点，其值应当以点和P、R点的插值给出代入第二、三类边界条件2021/6/1624图中O与R重合图中V与R点重合（2）第二类和第三类边界条件 2021/6/16252.4 差分方程差分方程对于具体地球物理问题的偏微分方程组，利用上述差分格式，可以给出偏导数的微商近似，进一步得到差分方程组。设f(P)是内部区域DI上定义的一个函数，设L(u)是一个微分算子，则以下表示了未知量u(P)的偏微分方程：边界条件表示为以下方程：设 和 分别表示区域D的北部节点和边界节点，则下式表示了以上偏微分方程的有限差分方程(finite-difference equations,or finite-difference scheme)：I()(),L u Pf PPDB()(),B u Pg PPDI()(),L Uf PPDB()(),B Ug PPD其中：以上差分方程是偏微分方程的有限差分近似，U是u的有限差分近似。差分方程要求U在所有节点上是u的很好的近似，并且方程所给出的有限差分近似解U是唯一的。IDBD2021/6/1626例子：对流方程（双曲型）的初值问题差分方程差分方程0,0(,0)(),uuaxR ttxu xg xxR(13)(14)假定以上问题的解u(x,t)是充分光滑的，由Taylor级数展开有：2021/6/1627利用（15）和（17）式，得到：如果u(x,t)是满足方程（13）的光滑解，则代入（20），可以看出，偏微分方程（13）在(xj,tn)处可以近似地用下面的方程来代替：其中是u(xj,tn)的近似值。（21）式称为逼近方程（13）的有限差分方程，简称差分方程。用到的节点如图所示，（21）式可以改写为便于计算的形式11(,)(,)(,)(,)()njnjnjnjnju x tu x tu xtu x tuuaaOhhtx(20)0njuuatx1100,1,2,n=0,1,2,nnnnjjjjuuuuahj(21)11()nnnnjjjjuuauu(22)2021/6/1628（22）式加上初始条件（14）的离散形式0,0,1,2,jjuj 11()nnnnjjjjuuauu(22)（23）就可以按时间逐层推进，算出各层的值。这里的“层”表示在直线t=n上网格点的整体。差分方程（22）和初始条件的离散形式（23）结合在一起构成了一个差分格式。（22）给出了根据初始条件（23）来确定(j=0,1,)的一个算法，因此有时也称差分方程（22）为一个差分格式。差分格式包含了初始条件、边界条件的离散。由第n个时间层推进到第n+1个时间层时，（22）提供了逐点直接计算n+1时的表达式，因此称（22）为显式格式，在计算第n+1层时只用到了n层的数据，前后仅联系到连个时间层，因此（21）式为两层格式，更明确地称其为两层显式格式。用（15）和（19），可以得到逼近方程（13）的另一差分格式11-1020,1,2,n=0,1,2,nnnnjjjjuuuuahj 11-1()2nnnnjjjjauuuu其中为网格比。此格式也是两层格式，称为中心差分格式。2021/6/1629用Taylor级数展开给出差分形式，用相应的差分形式逼近批未能微分方程（组），可以得到不同的差分格式，这一过程等价于用差商来近似微商得到相应的差分格式。不同的差分格式的精度和误差不同。2021/6/1630例子：牛顿冷却定律：温度高于周围环境的物体向周围媒质传递热量逐渐冷却时所遵循的规律。当物体表面与周围存在温度差时，单位时间从单位面积散失的热量与温度差成正比。Tair一阶常微分方程的数值解(,)dTf T tdt如果不能简单地通过积分求解，则需要利用数值方法求解。首先对时间和温度进行离散：0,()jjjttj t TT t 利用向前差分形式：1()jjjt tTTdTOtdtt得到以下的显式差分格式：1(,)jjjjTTtf T tTcap2021/6/1631利用该差分格式，我们可以计算任一时间和函数f的温度，但是随着计算的进行，误差O(dt)将会不断积累。1(,)jjjjTTtf T t 对于例子中的牛顿冷却问题，我们想知道液体的温度怎样随着时间和它与空气的温度差变化。设温度差为capairTTT为冷却的时间尺度，这时(,)/f T tT 常微分方程为dTTdt 初始条件：0,0TTt以上的差分格式为：1(1)jjjjttTTTT该方程的解析解为：0()exp(/)T tTt该有限差分格式的近似程度如何？怎样选择t才能得到稳定解？2021/6/1632差分方程1(1)jjtTT该有限差分格式是否在t0时收敛？它是否和解析解得性质一样？考虑初始条件0,0TTt得到10(1)tTT2210(1)=(1)ttTTT因此，对于第j个时刻0(1)jjtTT设t是到j时刻的总时间，/ttj 则：0(1)jjtTTj122011()1()12jjjjjjttTTjj 对上式利用二项式展开其中：!()!jjrjr r 我们希望知道当dt0时差分格式的结果，这时相当于j!(1)(2)(1)()!rjj jjjrjjr 2021/6/1633因此：!rjjrr 代入差分格式：代入差分格式：2120201,1!2!11,2!jjjtjtTTjjttTT 上式就是解析解0()exp(/)T tTt的Taylor展开。结论：对于牛顿冷却问题，当时间步长t趋于零时，差分格式给出的数值解收敛于解析给出的严格解。解析解是单调减小函数，数值解的性质怎样？保证数值解单调减小（Tj+1Tj）需要什么条件？122011()1()12jjjjjjttTTjj 得到得到!()!jjrjr r !(1)(2)(1)()!rjj jjjrjjr 2021/6/16341(1)jjtTT011t 0t 保证Tj+11e-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=4192021/6/165405101505100204060801005101512345678910112021/6/16552.5 应用实例应用实例南加州一次未来大地震的强地面运动的数值模拟2021/6/1656盆地效应2021/6/1657Cui,20132021/6/1658Cui,20132021/6/1659Cui,20132021/6/1660Cui,20132021/6/1661总结：总结：1、有限差分方法给出的数值解的精度取决于所用的差分形式（向前、向后、中心）。2、偏微分方程的显式有限差分格式通常是有条件稳定的，为了保证得到精确的数值解，最关键的是需要根据稳定性条件选取正确的空间和时间步长。若有不当之处，请指正，谢谢！若有不当之处，请指正，谢谢！