有限差分法

第四章有限差分方法4.1引言有限差分法：数值求解常微分方程或偏微分方程的方法。物理学和其他学科领域的许多问题在被分析研究之后,往往可以归结为常微分方程或偏微分方程的求解问题。一般说来，处理一个特定的物理问题，除了需要知道它满足的数学方程外，还应当同时知道这个问题的定解条件，然后才能设计出行之有效的计算方法来求解。有限差分法以变量离散取值后对应的函数值来近似微分方程中独立变量的连续取值。在有限差分方法中，我们放弃了微分方程中独立变量可以取连续值的特征，而关注独立变量离散取值后对应的函数值。但是从原则上说，这种方法仍然可以达到任意满意的计算精度。因为方程的连续数值解可以通过减小独立变量离散取值的间格，或者通过离散点上的函数值插值计算来近似得到。这种方法是随着计算机的诞生和应用而发展起来的。其计算格式和程序的设计都比较直观和简单，因而，它的实际应用已经构成了计算数学和计算物理的重要组成部分。1有限差分法的具体操作分为两个部分：(1)用差分代替微分方程中的微分，将连续变化的变量离散化，从而得到差分方程组的数学形式；(2)求解差分方程组。在第一步中，我们通过所谓的网络分割法，将函数定义域分成大量相邻而不重合的子区域。通常采用的是规则的分割方式。这样可以便于计算机自动实现和减少计算的复杂性。网络线划分的交点称为节点。若与某个节点P相邻的节点都是定义在场域内的节点，则P点称为正则节点；反之，若节点P有处在定义域外的相邻节点，则P点称为非正则节点。在第二步中，数值求解的关键就是要应用适当的计算方法，求得特定问题在所有这些节点上的离散近似值。有限差分法的差分格式：一个函数在x点上的一阶和二阶微商，可以近似地用它所临近的两点上的函数值的差分来表示。如对一个单变量函数f(x)，x为定义在区间[a,b]的连续变量。以步长h=Δx将[a,b]区间离散化，我们得到一系列节点x=a,x=x+h,x=x+h=a+212132Δx,...,x=x+h=b,然后求出f(x)在这些点上的近似值。显然步长h越小，近似解的精度就越好。与节点相邻的节点有和n+1nixhxi−hxi+，因此在点可以构造如下形式的差值：ix),()(iixfhxf−+节点的一阶向前差分xi),()(hxfxfii−−节点的一阶向后差分xi).()(hxfhxfii−−+节点的一阶中心差分xi2与点相邻两点的泰勒展开式可以写为xifxhfxhfxhfxhfxhfxiiiiii()()()()!()!()......''''''''''−=−+−+−234234.(4.1.1)fxhfxhfxhfxhfxhfxiiiiii()()()()!()!()......''''''''''+=+++++234234.(4.1.2)（4.1.1）-(4.1.2)，并忽略h的平方和更高阶的项得到一阶微分的中心差商表示：fxfxhfxhhiii'()()()≈+−−2.(4.1.3)利用(4.1.1)和(4.1.2)式我们还可以得到一阶微分的向前，向后一阶差商表示fxfxhfxhiii'()()()≈+−,(4.1.4)fxfxfxhhiii'()()()≈−−.(4.1.5)将(4.1.1)和(4.1.2)式相加，忽略h的立方及更高阶的项得到二阶微分的中心差商表示fxfxhfxfxhhiiii''()()()()≈+−+−22.(4.1.6)3利用(4.1.3)~(4.1.6)式，我们就可以构造出微分方程的差分格式。这里要指出的是：在构造差分格式时，究竟应该选择向前，向后还是中间差分或差商来代替微分方程中的微分或微商，应当根据由此得到的差分方程解的稳定性和收敛性来考虑。同时兼顾到差分格式的简单和求解的方便。上述差分步骤应用于偏微分：例如，对于的情况，拉普拉斯算符在0点作用在此函数上的值ffxy=(,)22222fffxy∂∂∂∂⎛⎞⎛⎞⎟，也可以用临近的点上的函数值来表示出来。(见图4.1.1,且∇=+⎜⎟⎜⎝⎠⎝⎠hhhhh====4321时)).(!424444422043212yfxfhhffffff∂∂∂∂+−−+++≈∇(4.1.7)j-1ji+1j-1(i-1,j-1)(i,j-1)h1h34h4(i+1,j-1)(i+1,j)(i,j)j(i-1,j)301h22(i+1,j+1)(i,j+1)j+1(i-1,j+1)图4.1.1节点0及邻近节点4对微分方程数值求解的误差的来源：（1）方法误差（或截断误差）。这是由于采用的计算方法所引起的误差。例如上面我们介绍的差商表示中，采用的泰勒展开式展开到第1+n项时的截断误差阶数为()1+nhO。具体方法的误差阶数取决于在离散化时的近似阶数。因此若改进算法就可以减小截断误差。（2）舍入误差（或计算误差）。这是由于计算机的有限字长而造成数据在计算机中的表示出现误差。在计算机运算的过程中，随着运算次数的增加舍入误差会积累得很大。如果在多次运算后，舍入误差的精度影响是有限的，那么这个算法是稳定的，否则是不稳定的。不稳定的算法是不能用的。本书中我们将略去对差分法稳定性和收敛性理论的讨论，尽管这方面的内容是相当重要的。以下的讨论中所讲到的各种差分格式，我们均假定求解方法满足稳定性和收敛性的要求。54.2有限差分法和偏微分方程利用上节所介绍的微分的差分表示，我们就很容易地将微分方程离散化为差分方程组的形式。但是由差分方程所得的解完全取决于待求微分方程的特性。正如我们在物理上所知道的，边界条件的情况变化将会引起差分方程组的不同。在求解微分方程中，我们会遇到两类问题：一类是初始值问题；另一类是边值条件的问题。在初始值问题中，部分边界上的函数值和部分的函数偏导值是给定的。通常在这类问题中的独立变量之一是时间t。在边界值问题中，边界上的信息是给定的。本书中我们仅讨论后一类问题。假定某方程形式上可以写为：Lqφ=.(4.2.1)其中L为含偏微商的算符．它的边界条件一般可写为：).(|)(|21sgnsgGG=+∂∂φφ(4.2.2)G表示场域D的边界，为边界上s点的逐点函数。)(),(21sgsg6三类边界条件：（1）第一类边界条件，或称为狄利克莱(Dirichlet)问题(,)gg1200=≠。).(|sgG=φ(4.2.3)（2）第二类边界条件，或称诺伊曼(Neumann)问题(,)gg1200≠=。).(|sgnG=∂∂φ(4.2.4)（3）第三类边界条件，或称混合问题(,)gg1200≠≠。对于算符L为斯杜-刘维尔(Sturm-Liouville)算符的特定情况，即Lpf≡−∇∇+().(4.2.5)公式中的p和f是给定的函数。我们将会得到一类很重要的微分方程。它是在流体力学、等离子物理、天体物理…等学科中，势函数起关键作用的许多问题当中的基本方程。当p=1,f=0时，我们得到（4.2.1）式的特殊情况—泊松(Poisson)方程。我们现在考虑方程(4.2.1)中p为常数的二维的情况，我们可以得到下面的方程：∂φ∂∂φ∂φ2222xyfxyqxy++=(,)(,).(4.2.6)设函数φ在区域D内满足方程（4.2.6）式(区域的边界为G)。D7区域的离散化：D即通过任意的网络划分方法把区域离散为许许多多的小单元。原则上讲这种网格分割是可以任意的，但是在实际应用中，常常是根据边界G的形状，采用昀简单，昀有规律，和边界的拟合程度昀佳的方法来分割。常用的有正方形分割法和矩形分割法（如图4.2.1）。有时也用三角形分割法（见图4.2.2）。对圆形区域，应用图（4.2.3）所示的极网络格式也许更方便些。这些网络单元通常称为元素，网络点称为节点。DD0xyD0xy4.2.1求解区域的矩形分割。4.2.2求解区域的三角形分割。4.2.3求解区域的极网络分割。8用节点上的函数值来表示节点上偏导的数值。设区域内部某节点0附近的各节点如图4.1.1所示。这里我们取步长h不相等的昀一般情况。以φφφφφ01234分别代表在节点0，1，2，3，4处,,,,φ的函数值。如前所述，0点的一阶偏导数可以通过先前或向后的差商，由1和3节点近似写出∂φ∂φφxh⎛⎝⎜⎞⎠⎟≈−0101.（4.2.7）或∂φ∂φφxh⎛⎝⎜⎞⎠⎟≈−0033.（4.2.8）显然这种单侧差商的误差较大。如果要寻求更精确的差分格式，我们可以引入待定常数αβ,，由φ1和φ3的泰勒展开，构造出如下的关系式：....)(21)()()(23210223100301++⎟⎟⎠⎞⎜⎜⎝⎛+−⎟⎠⎞⎜⎝⎛=−+−hhxhhxβα∂φ∂βα∂∂φφφβφφα（4.2.9）令∂φ∂220x⎛⎝⎜⎞⎠⎟项的系数为零，则得到α和β之间应当满足αβ=−hh3212.（4.2.10）9将公式（4.2.10）带入（4.2.9），并舍去高阶项，得到∂φ∂x⎛⎝⎜⎞⎠⎟0的另一个差分表达式∂φ∂φφφφxhhhhhh⎛⎝⎜⎞⎠⎟≈−−−+0321012301313()()().（4.2.11）当选用等步距时，上式成为hhhx13==∂φ∂φφxhx⎛⎝⎜⎞⎠⎟≈−0132.（4.2.12）（此为中心差商表达式。）二阶偏导数的差分表达式（“五点格式”或“菱形格式”）在（4.2.9）式中，如果令∂φ∂x⎛⎝⎜⎞⎠⎟0的系数为零，则有α和β间存在关系式：αβ=hh31.（4.2.13）将上式（4.2.13）代入（4.2.9）式中，并忽略h三阶以上的高次项，则得到表达式：∂φ∂φφφφ22031013013132xhhhhhh⎛⎝⎜⎞⎠⎟≈−+−+()()().（4.2.14）当用等步长hhhx13==时，上式成为10∂φ∂φφφ220103（4.2.15）22xhx⎛⎝⎜⎞⎠⎟≈−+.用完全相同的计算方法，我们可以推导出∂φ∂220y⎛⎝⎜⎞⎠⎟的差分表达式：∂φ∂φφφφ22042024024242yhhhhhh⎛⎝⎜⎞⎠⎟≈−+−+()()().（4.2.16）当采用等步长时，有hhhy24==∂φ∂φφφ22020422yhy⎛⎝⎜⎞⎠⎟≈−+.（4.2.17）将公式（4.2.14）和（4.2.16）两式代入方程（4.2.13），我们就得到该方程的差分表达式为0004242042024313103101302)()()()()()(2)(qfhhhhhhhhhhhh=+⎥⎦⎤⎢⎣⎡+−+−++−+−=∇φφφφφφφφφφ（4.2.18）如果在x和y方向的步长分别相等，即hhhx13==和hhhy24=时，则上式化为=φφφφφφφ1032204200022−++−++=hhfqxy，（4.2.19）一般可以用角标来表示节点的标记，将上式写为111212211211hhfqxijijijyijijijijijij()(),,,,,,,,,，(4.2.20）φφφφφφφ+−+−−++−++=特别是当hhhxy==时，我们得到：,,（4.2.21）φφφφφijijijijijijijhfhq+−+−++++−=1111224,,,,,,()f=0的时候方程（4.2.6）为泊松方程，由（4.2.21）式得到,，（4.2.22）φφφφφijijijijijijhq+−+−+++−=111124,,,,,fq==0的时为拉普拉斯方程，从（4.2.21）式得φφφφφijijijijij+−+−+++−=111140,,,,,，（4.2.23）边界条件的离散化的处理：若场域的网络节点都落在边界G上，则显然无需再做处理。但是在一般情况下，边界G是不规则的。网络节点不可能全部都落在边界G上。对(4.2.3)式给出的第一类边界条件，通常有两种处理办法。一种是所谓的直接转移法，如果0节点靠近边界，则取昀靠近0点的边界节点上的函数作为0点的函数值。这是一种比较粗糙的近似。另一种方法是较为精确的线性插值法。对第二、三类边界条件也可以用插值法求出临近边界节点上的函数值。12（1）第一类边界条件在二维情况下，第一类边界条件如公式(4.2.3)所示。如果网格的边界节点恰好落在边界G上，则显然无需再做近似处理，边界节点的函数值就等于边界条件(4.2.3)给