第4章有限差分法目录第4章有限差分法14.1有限差分法的原理24.1.1定解区域离散化34.1.2微分方程离散化34.1.3初边值条件离散化.44.1.4解差分方程组.44.2椭圆型方程的有限差分法54.2.1微分方程的离散化61、五点差分网格62、九点差分网格.74.2.2边值条件的离散化81、第一边界条件离散化82、第三边界条件离散化94.3抛物型方程的有限差分法104.3.1 维抛物型方程的差分格式101、向前差分格式112、向后差分格式113、六点对称格式124、Richardson 格式124.3.2高维方程的差分格式131、二维问题差分方法132、分裂格式算法144.4双曲型方程的有限差分法164.4.1 一阶线性双曲型方程的差分格式171、一阶线性双曲型方程初值问题172、一阶线性双曲型方程初边值问题184.4.2二阶线性双曲型方程的差分格式191、波动方程及其特征192、二阶方程的差分格式203、定解条件的处理21参考文献23第4章有限差分法刚刚介绍的变分法和加权余量法，都是将实际问题的控制方程转化为积分形式进行求解， 而有限差分方法不同，它是用网格覆盖求解区域和时间域，用差分近似替代控制方程中的微 分，把求解微分方程的问题改换成求解代数方程问题，进行近似求解。有限差分法是一种经典数值方法，当网格细化或节点较多时，近似解的精度可以得到改 进，对于几何形状规则且边界条件简单的问题，差分法有很高的求解精度和收敛性，但用于 求解几何形状和边界条件复杂的问题时，它的精度将会降低，甚至产生困难。有限差分法在 流体力学和爆炸力学中得到广泛应用，也曾用以求解弹性力学边值问题。目前，弹性静力学 一般都不用差分法进行求解了，但是对于弹性动力学等初值问题，仍然要借助于有限差分法 予以解决。有限差分法的关键是构造差分格式，保证解的收敛性和解的稳定性。这里简要介绍椭圆 型方程、抛物型方程和双曲型方程的差分格式构造方法。4.1有限差分法的原理差分法是微分方程的一种近似数值解法。它不是去寻求函数式的解答，而是去求出函数 在一些网格结点上的数值。具体地讲，差分法就是先将求解域离散化，在这些离散点上用差 商近似地代替导数，从而把基本方程和边界条件近似地改用差分方程来表示，把求解微分方 程的问题改换为求解代数方程的问题。考虑一维热传导方程的第一类初边值问题Qu =aQtQ 2UQx2+ f(x) u(x,0)=甲(x) u(0,t) = u(1,t) = 0(x,t) e Qx e (0,1)(4-1) t e 0,T其中，。=(x,t)| x e 0,1,t e 0,T，a 为正常数。假定f(x)和中(x)在相应区域光滑，并且在x = 0,l满足相容条件，使上述问题有唯 一充分光滑的解。用差分法求解上述问题时，通常情况下，首先采用矩形网格对求解区域进行离散化，接 着，利用差商或Taylor展式对微分方程和初边值条件等连续性方程，在网格节点处做离散 化处理，最后解差分方程组，获得所有网格节点处的解的近似值，用近似研究微分方程初边 值问题的解。4.1.1定解区域离散化由于有限差分法求解偏微分方程问题必须把连续问题进行离散化，因此首先要对求解区域进行网格剖分。对于不同的求解问题，求解区域也各不相同。对于问题（4-1）而言，如 4-1所示，设N是一给定的正整数，取空间步长Ax =i/n，再取一定的时间步长At =t/M，做图4-1差分法矩形网格剖分4.1.2微分方程离散化由微分学知道，函数的导数是函数的增量与自变量增量之比的极限，即u，- lim u(x + Ax) - u(x) _ jm u(x) 一 u(x - Ax)Axr 0Axaxt 0Au- lim u(x +Ax)- u(x)Ax t 0Ai .1u(x + Ax) - u(x)u(x) - u(x - Ax)-limAx t 0 Ax LAxAx_. u(x + Ax) - 2u(x) + u(x - Ax)-lim Ax t0(Ax)2当Ax很小时，u，可以近似地用差商u(x + Ax) - u(x)Axu(x) - u(x - Ax) 或Ax代替，u 可以近似地用二阶差商u(x + Ax) - 2u(x) + u(x - Ax)（Ax）2代替，从而一个微分方程可近似地用一个差分方程来代替。因此，对于问题（4-1）而言， 其一维热传导方程du _ 2 d2U就dx2可以用方程u（x，t +At）一u（x,t）u（x+Ax,t） 一2u（x,t） +u（x-Ax,t） _ a2 （4-2）At（Ax）2近似代替，其截断误差为O（At）+ O（Ax）2）。差分格式（4-2）称为解热传导方程（4-1）的古典显格式，它是一个两层4点格式，所用到的节点图式如图4-2所示。（j，k+i）（j-1,k）（j,k）（j+1,k）图4-2古典显式差分格式节点图在有限差分法中，为简单起见，通常将节点（x,t）处的函数值U（X,t）简记为uk，依 j kj次，节点（X ,t）处的函数值u（x +Ax,t）简记为uk，节点（X ,t ）处的函数值j+1 kjj+1j k+1u（xj,t +At）简记为uk+1。这样的话，式（4-2）就可以写成Uk +1 uk uk 2uk + uk-j_ a2jj-1 + f k (4-3)Th2j4.1.3初边值条件离散化对于问题（4-1）中的初边值条件，在网格节点处取值，有u0 _ q)(x ), uk _ uk _ 0, jj0 Nuk _ uk _ 0, k _ 1,2, M.0N4.1.4解差分方程组对于问题（4-1），将微分方程和初边值条件的离散化方程联立，得差分方程组Uk +1=uk + ruk2uk + uk +T f kjjj 1jj +1j=ruk + (12 r)uk+ ruk+ T f k,j 1jj +1j（4-4）U0 =q)(x),uk=uk =0,jj0N、r = a2T / h2,j = 1,2,N 1, k = 1,2,M 1这样，利用k=0层上的数值，即可逐点算出k=1层上的全部离散点处解的近似值，再利 用k层上的节点值算出k+1层上的节点值。这样逐层逐点地计算数值解是非常方便的。同样地，若给定的是初值问题，仍采用前面的古典显格式和初值离散方法，也容易写出 其差分方程组，但其计算求解范围与初边值问题情形有所不同。此时，利用k=0时间层上的 节点近似值只能计算出k=1时间层上关于j=1,2, N-1处的节点值近似值，依次下去，纯初值 问题的差分解被局限于图4-1中0A1上的网格节点处。4.2椭圆型方程的有限差分法这里以Possion方程边值问题为例，讨论椭圆型方程的差分格式建立，边界条件处理，以及 差分格式的收敛性问题的讨论。考虑Possion方程仞U仞UAu =+ =f（x，y） （x，y） e D （4-5） dx2Qy 2其中D是x-y平面内一有界区域，其边界用r表示，假定它是分段光滑的曲线所组成，如图 4-3所示，对求解域D剖分矩形网格，在平面上取一固定点（xo,y0），设x轴方向的步长为h， y轴方向的步长为t，过该点做垂直于x轴和y轴的平行线x = x + ih,i = 0,1,N ，y = y0 + jT,j = 0,1 ,N2，称网格线的交点（七，）为网格点。位于d内部的网格点称为内部节点，简称内点，它们的全体记为Qh ；网格线与区域边界r的交点称为边界节点，简称边界点，它们的全体记为卜。为了便于后面讨论，将内点分为两类。若某内点的网格线上四个相邻的网格点都是内点， 则称它为正则内点，否则称为非正则内点。4.2.1微分方程的离散化1、五点差分网格现在我们来构造Possion方程的五点差分格式。为简单起见，把点（xi,yj）记为（i,j）。参见图4-4所示，我们这里采用Taylor展开方法给出差分格式。设（i,j）为内节点，则有_1lh21()d2ujh2 d4ujh4 6ujZ1 、uj 2uj + uj = i + 1 i + 1i + 0(h6) (4-6)i+1ii-18x212 8x4360 8x61_1l h21 2uj + uj-1)=划 + 竖竺 + 刘竺 + 0(h6) (4-7) ii i8y212 8y4360 8y622d 2u 利用式（4-6）和（4-7）中沿X 和 y方向的二阶中心差商直接代替方程（4-5）中的应合2U和，就得到（4-8）Cjj 2uj + uj )+4 Cjj+i 2uj + uj-i) = fj h2 i+1 ii-1 h2 ii ii12由于在建立差分方程（4-5）时用到了节点（i，j）及其相邻的四个节点，所以称式（4-8） 为五点差分格式。差分方程相当于在节点（i，j）分别沿x、y方向的二阶中心差商代替二阶 导数的结果。由式（4-6）和（4-7）可知五点差分格式的截断误差为Rj = Auj A uj = 0（h2） ii h i其中，h=max（hi，h2）o截断误差反映了差分算子对微分算子的相容逼近，也反映了差 分解的局部误差状况。（i,j+1）II（i-1,j） II（i,j）（i+1,j）（i,jT）图4-4五点差分格式节点图图4-3求解域离散化如果h1h2，即用正方形网格，则差分方程化为A hUj1 (h2Tuj + ujii-1+ u ji+1+ uj+1i+ u j-1)=i(4-9)形式特别简单，若f=0,则u.4 i-1,j+ui + 1,j+ui,j +1+ ui,j T(4-10)差分解在节点（i，j）的值等于其周围四点值的平均。2、九点差分网格为了提高差分格式的精度，考虑九点差分格式，为推导简单起见，令h = t， 编号如图4-5所示。这里我们按Taylor展开推导差分格式，得j+1 +uj-1 +uj+1 +uj-1 -4u j）节点及 uj三i2h2i+1i+1i-1i-1 i.h2 （ 54uA uj+i+12 5x46 5 4u +5x2y25y4 ，j)(4-11)h4+ 一155x4 Qy2155 6u56u+5x25y45y6 J 、4, j)+ 0(h6)uj =ij +uj +uj+1+uj-1-4uj) h2i+1i-1 i i i(54u+ 12 5x45y4.h254u