有限差分法基础课件

第二章有限差分法主讲人：胡才博中国科学院大学地球科学学院中国科学院计算地球动力学重点实验室第二章有限差分法主讲人：胡才博第二章有限差分法2.1有限差分法基础2.2网格剖分2.3差分格式2.4差分方程2.5应用实例第二章有限差分法2.1有限差分法基础1.地球内部介质，不仅存在纵向非均匀结构（一维地球模型），也存在横向非均匀结构（不同块体、断层系统）；2.几何模型也呈现出相当的复杂性；3.另外，边界条件和初始条件对于不同问题具有特殊性。解析方法的局限性Hu,C.,Y.Cai,andZ.Wang(2012),Effectsoflargehistoricalearthquakes,viscousrelaxation,andtectonicloadingonthe2008Wenchuanearthquake,JournalofGeophysicalResearch,117,B06410,doi:10.1029/2011JB009046.(SCI,IF:3.303)汶川大地震的动力学成因1.地球内部介质，不仅存在纵向非均匀结构（一维地球模型），对于存在复杂介质和几何、特殊边界条件和初始条件的实际地质问题，一般不存在解析解，需要近似的数值求解方法。有限差分方法是地球物理方法中最常见的一种。有限差分方法（FiniteDifferenceMethod,FDM）是计算机数值模拟最早采用的方法，至今仍被广泛使用。有限差分方法的基本特点该方法是一种直接将微分问题变为代数问题的近似数值解法，数学概念直观，表达简单，是发展较早且比较成熟的数值方法。2.1有限差分法基础对于存在复杂介质和几何、特殊边界条件和初始条件的实际地质问题有限差分法以变量离散取值后对应的函数值来近似微分方程中独立变量的连续取值。我们放弃了微分方程中独立变量可以取连续值的特征，而关注独立变量离散取值后对应的函数值。有限差分法的具体操作分为两个部分：(1)用差分代替微分方程中的微分，将连续变化的变量离散化，从而得到差分方程组的数学形式；(2)求解差分方程组。有限差分方法的基本原理该方法将求解域划分为差分网格，用有限个网格节点代替连续的求解域。有限差分方法以Taylor级数展开等方法，把控制方程中的导数用网格节点上的函数值的差商代替进行离散，从而建立以网格节点上的函数值为未知数的代数方程组。2.1有限差分法基础有限差分方法的基本原理该方法将求解域划分为差分网格，用有限个有限差分法的主要内容建立地球物理问题的离散有限差分模型（1）如何根据问题的特点将定解区域做网格划分；（2）如何在所有网格节点上用有限差分格式对导数求近似，对函数、初始条件和边界条件求近似；（3）如何把原方程离散化为代数方程组，即有限差分方程组。2.从理论上研究有限差分模型的形态，以保证计算过程的可行性和计算结果的正确性（1）解的相容性；（2）解的稳定性；（3）解的收敛性。3.如何数值求解差分方程组有限差分法的主要内容建立地球物理问题的离散有限差分模型2.从网格剖分就是研究区域和边界的离散化1.矩形分割2.三角形分割3.极网格分割2.2网格剖分网格剖分就是研究区域和边界的离散化2.2网格剖分对地球物理问题的连续求解区域通过网格划分离散为空间上得一系列网格点，接下来需要利用一定的差分格式对偏微分方程组中的导数用差商进行近似，从而将偏微分方程组离散化为差分方程组。对于函数f(x)，通常意义下的导数（微商）定义为：当dx→0时，以上三种形式都是微商的正确定义。如果dx是有限的，如何给出微商的近似定义？对地球物理问题的连续求解区域通过网格划分离散为空间上得一系列2.3差分格式用Taylor级数展开可以给出微商的近似形式。对于连续函数f(x)，它在相邻点上的值f(x+Δx)和f(x-Δx)可以用Taylor级数展开为dx变为Δx如果Δx很小，f(x)可微，则以上级数收敛。次数越高，收敛级数的项的绝对值越小。由（1）得到，（1）（2）（3）（4）2.3差分格式用Taylor级数展开可以给出微商的近似形式式中的O(x)项表示忽略掉的所有项中的最大项的量级是Δx，也就是说，忽略掉这些项带来的误差中的最大项和Δx成正比。由（4）给出导数的一阶精度（firstorderaccurate）近似为：（4）（5）（5）式称为向前差分格式（forward-differenceformula）由（2）式得到由（7）式得到导数的另一个一阶精度近似：（6）（7）（8）（8）式称为向后差分形式（backward-differenceformula）。式中的O(x)项表示忽略掉的所有项中的最大项的量级是Δx，也（1）式减去（2）式，得到：（9）式中的O(Δx2)项表示忽略掉这些项带来的误差中的最大项和Δx2成正比。由（9）式得到导数的二阶精度（secondorderaccurate）近似为：（10）式称为中心差分形式（central-differenceformula）。（9）（10）（1）（2）（1）式减去（2）式，得到：（9）式中的O(Δx2)项表示忽（1）（2）（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）（1）式和（2）式相加，得到：（12）式称为二阶总结：1、向前差分形式：2、向后差分形式：3、中心差分形式：单侧，一阶精度单侧，一阶精度对称，二阶精度对于二阶导数二阶精度对一阶导数总结：1、向前差分形式：2、向后差分形式：3、中心差分形式：定解问题的有限差分解法1．离散x=ih,y=jh,i=0,±1,±2,….±n,h:步长（正方形的边长）2．根据泰勒级数建立差商格式：对于一维情况:在x处的一阶导数可以用3.建立和求解差商方程组3.建立和求解差商方程组差分格式的另一种推导为了寻求更精确的差分格式，我们引入两个待定常数，由泰勒展开，构造如下关系式(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差分格式的另一种推导为了寻求更精确的差分格式，我们引入两个待为了寻求更精确的差分格式，我们引入两个待定常数，由泰勒展开，构造如下关系式回代（1）中，舍去高阶项一阶偏导数中心差分的推导（1）为了寻求更精确的差分格式，我们引入两个待定常数，由泰勒展开，（1）二阶偏导数差分的推导回代（1）中，舍去高阶项二阶偏导数差分公式（1）二阶偏导数差分的推导回代（1）中，舍去高阶项二阶偏导数一个例子：等步长一个例子：等步长一个例子：局部节点离散化方程总体节点离散化方程总体节点离散化方程f=0时，变为泊松方程f=q=0时，变为拉普拉斯方程一个例子：局部节点离散化方程总体节点离散化方程总体节点离散化（1）第一类边界条件（a）直接转移法在图中网格是按正方形分割，步长为h。0点为靠近边界G的一个网格节点，1和2为边界节点。我们取最靠近0点的边界节点1上的函数值作为0点的函数值。即取φ0≈φ1。这种方法称为直接转移法，又称为零次插值法。

边界条件的离散化的处理狄利克莱问题（1）第一类边界条件（a）直接转移法边界条件的离散化的处（b）线性插值法先判断x方向的边界节点1和y方向的边界节点2哪一个更靠近0点。如果1更靠近0点，则可以用x方向的线性插值给出0点的函数值如果2更靠近0点，则可以用x方向的线性插值给出0点的函数值（b）线性插值法如果1更靠近0点，则可以用x方向的线性插值（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变步长二次偏导数（c）双向插值法(i,j)(i+1,j)(i-1,j)01（2）第二类和第三类边界条件对于点O过O点向边界G做垂线PQ交边界于Q，交网线段VR于P，OP=ah，PR=bh，VP=ch因为P一般不是节点，其值应当以点和P、R点的插值给出代入第二、三类边界条件（2）第二类和第三类边界条件对于点O过O点向边界G做垂线图中O与R重合图中V与R点重合（2）第二类和第三类边界条件图中O与R重合图中V与R点重合（2）第二类和第三类边界条件2.4差分方程对于具体地球物理问题的偏微分方程组，利用上述差分格式，可以给出偏导数的微商近似，进一步得到差分方程组。设f(P)是内部区域DI上定义的一个函数，设L(u)是一个微分算子，则以下表示了未知量u(P)的偏微分方程：边界条件表示为以下方程：设和分别表示区域D的北部节点和边界节点，则下式表示了以上偏微分方程的有限差分方程(finite-differenceequations,orfinite-differencescheme)：其中：以上差分方程是偏微分方程的有限差分近似，U是u的有限差分近似。差分方程要求U在所有节点上是u的很好的近似，并且方程所给出的有限差分近似解U是唯一的。2.4差分方程对于具体地球物理问题的偏微分方程组，利用上述例子：对流方程（双曲型）的初值问题差分方程(13)(14)假定以上问题的解u(x,t)是充分光滑的，由Taylor级数展开有：例子：对流方程（双曲型）的初值问题差分方程(13)(14)假利用（15）和（17）式，得到：如果u(x,t)是满足方程（13）的光滑解，则代入（20），可以看出，偏微分方程（13）在(xj,tn)处可以近似地用下面的方程来代替：其中是u(xj,tn)的近似值。（21）式称为逼近方程（13）的有限差分方程，简称差分方程。用到的节点如图所示，（21）式可以改写为便于计算的形式(20)(21)(22)利用（15）和（17）式，得到：如果u(x,t)是满足方程（（22）式加上初始条件（14）的离散形式(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）的另一差分格式其中λ为网格比。此格式也是两层格式，称为中心差分格式。（22）式加上初始条件（14）的离散形式(22)（23）就可用Taylor级数展开给出差分形式，用相应的差分形式逼近批未能微分方程（组），可以得到不同的差分格式，这一过程等价于用差商来近似微商得到相应的差分格式。不同的差分格式的精度和误差不同。用Taylor级数展开给出差分形式，用相应的差分形式逼近批未例子：牛顿冷却定律：温度高于周围环境的物体向周围媒质传递热量逐渐冷却时所遵循的规律。当物体表面与周围存在温度差时，单位时间从单位面积散失的热量与温度差成正比。Tair一阶常微分方程的数值解如果不能简单地通过积分求解，则需要利用数值方法求解。首先对时间和温度进行离散：利用向前差分形式：得到以下的显式差分格式：Tcap例子：牛顿冷却定律：温度高于周围环境的物体向周围媒质传递热量利用该差分格式，我们可以计算任一时间和函数f的温度，但是随着计算的进行，误差O(dt)将会不断积累。对于例子中的牛顿冷却问题，我们想知道液体的温度怎样随着时间和它与空气的温度差变化。设温度差为τ为冷却的时间尺度，这时常微分方程为初始条件：以上的差分格式为：该方程的解析解为：该有限差分格式的近似程度如何？怎样选择Δt才能得到稳定解？利用该差分格式，我们可以计算任一时间和函数f的温度，但是随着差分方程该有限差分格式是否在Δt→0时收敛？它是否和解析解得性质一样？考虑初始条件得到因此，对于第j个时刻设t是到j时刻的总时间，则：对上式利用二项式展开其中：我们希望知道当dt→0时差分格式的结果，这时相当于j→∞差分方程该有限差分格式是否在Δt→0时收敛？它是否和解析解得因此：代入差分格式：上式就是解析解的Taylor展开。结论：对于牛顿冷却问题，当时间步长Δt趋于零时，差分格式给出的数值解收敛于解析给出的严格解。解析解是单调减小函数，数值解的性质怎样？保证数值解单调减小（Tj+1<Tj）需要什么条件？得到因此：代入差分格式：上式就是解析解的Taylor展开。结论：保证Tj+1<Tj的条件是：因此，数值解只对一定取值范围的dt是单调减小的。如果数值解振荡，但是结果会收敛如果数值解振荡并且发散。因此，数值结果是有条件稳定的。保证Tj+1<Tj的条件是：因此，数值解只对一定取值范围的d差分方程fori=1:nt+1xt(i)=(i-1)*dt;T1(i)=t0*exp(-xt(i)/tau);endplot(xt,T1,'r.-');holdon

set(gca,'DataAspectRatio',[(max(xt)-min(xt))/(max(T)-min(T))/311]);xlabel('Time(s)','Fontname','timesnewroman','FontSize',14);

ylabel('Temperature','Fontname','timesnewroman','FontSize',14);%Malab-1Dclear;clc;figure('color','w');nt=8;%totaltimestepst0=1;%initialtemperaturetau=0.7;%timeconstantdt=1.25;%timeintervalT(1)=t0;fori=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.-');holdon解析解差分解差分方程fori=1:nt+1%Malab-1D解析解差结果振荡但是收敛，前期结果不准确结果收敛，前期结果不准确结果振荡但是收敛，结果收敛，结果收敛，误差逐渐缩小结果收敛，结果基本吻合解析解结果收敛，结果收敛，结果振荡，不收敛，计算结果错误结果振荡，不收敛，计算结果错误差分格式的性质分析差分格式的性质分析有限差分法基础课件有限差分法基础课件差分格式的性质分析差分格式的性质分析差分格式的性质分析稳定性（stability）：如果偏微分方程的严格解析解有界，差分格式给出的解也有界，称该差分格式是稳定的；如果差分格式给出的解是无界的，则称该差分格式是不稳定的。如果差分格式给出的解对于所有的时间步长和空间步长取值都是有界的，则称该差分格式是无条件稳定的；如果只是对时间步长和空间步长的部分取值有界，称它是有条件稳定的；如果对于所有的时间步长和空间步长取值都是无界的，则称差分格式是无条件非稳定的。稳定性反映了差分格式在计算中控制误差传递的能力，对偏微分方程有限差分方法研究具有重要意义。例子：对流方程（双曲型）的初值问题差分格式的性质分析稳定性（stability）：如果偏微分方有限差分法基础课件收敛性（convergence）：如果当时间和空间步长趋于零时，FDE解趋于PDE解，称该差分格式是收敛的。如果则称该差分格式是收敛的。对连续形式的偏微分方程进行有限差分离散时，差分格式对最终计算结果有重要影响。收敛性表示当差分网格无限细化时，差分方程的解是否具有无限逼近偏微分方程的解的能力。收敛性的判断决定了一个差分格式能否被用来离散偏微分方程，不收敛的差分格式无疑是无法得到偏微分方程的真实解的。收敛性（convergence）：如果当时间和空间步长趋于零利用向前差分，得到以下差分格式：当时间步长dt趋近于零时，差分格式的近似解趋近于方程的解析解说明：该差分格式是收敛的。收敛性实例利用向前差分，得到以下差分格式：当时间步长dt趋近于零时，差很重要的一点：相容性是差分格式的性质，它将差分格式和原始的偏微分方程联系在一起。稳定性和收敛性则是差分格式给出的数值解的性质。一般来讲，分析稳定性和相容性比较容易，证明收敛性有时是很困难的数学问题。Lax等价定理（Laxequivalencetheorem）：如果逼近一个给定问题的差分格式是相容的，那么该差分格式的收敛性与稳定性互为充分必要条件。该定理表明，对于一个具有相容性的差分格式，它的收敛性与稳定性是等价的。如果格式不稳定，则也不收敛。由于判断差分格式的收敛性相对比较困难，该定理提供了通过判断差分格式稳定性来确定收敛性的方法。因此，一般不再讨论收敛性问题，而是讨论差分格式的稳定性。该定理将收敛性与相容性和稳定性联系在一起，是很有用的一个定理。差分格式的性质分析很重要的一点：相容性是差分格式的性质，它将差分格式和原始的偏拉普拉斯方程U=0U=0U=100U=0拉普拉斯方程的有限差分法问题：2.5应用实例拉普拉斯方程U=0U=0U=100U=0拉普拉斯方程的问题：U=0U=0U=100U=0214365871091211141315U=0U=0U=100U=02143658710912111组建A和B矩阵，求解线性方程组得到X组建A和B矩阵，求解线性方程组得到X%Matlab2Dclear;clc;figure('color','w');

a=zeros(135,135);fori=1:135a(i,i)=1;end;fori=1:7a(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;endfori=1:7a(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;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;fori=2:14a(i,i-1)=-0.25;a(i,i+1)=-0.25;a(i,i+15)=-0.25;endfori=122:134a(i,i-1)=-0.25;a(i,i+1)=-0.25;a(i,i-15)=-0.25;endfori=1:7forj=2:14;a(15*i+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;endendb=a^(-1);