有限差分法

有限差分方法许多物理现象或其运动、演化过程可以用一个微分方程的定解问题来描述，例如无限长细弦的自由振动向题可归结成二阶双曲型方程的初值问题，而弦对平衡位置的偏移就是方程的解。但是，绝大多数偏微分方程定解问题的解通常并不能用显式的公式来表达，有时即使可用公式表示，也往往因为过于复杂，而需要采用各种近似方法来计算它的解。差分方法就是求解（偏）微分方程定解问题的常用近似方法之一。Courant,Friedrichs,Lewy(1928)首次对偏微分方程的差分方法作了完整的论述。第二次世界大战以来，快速电子计算机的诞生为差分方法提供了强有力的工具，从而促使这一数值分析方法迅速地发展起来。有限差分法的基本思想是用离散的、只含有限个未知数的差分方程去代替连续变量的微分方程和定解条件。对于求解的偏微分方程定解问题，有限差分方法的主要步骤如下：利用网格线将定解区域化为离散点集；在此基础上，通过适当的途径将微分方程离散化为差分方程，并将定解条件离散化，这一过程叫做构造差分格式，不同的离散化途径一般会得到不同的差分格式；建立差分格式后，就把原来的偏微分方程定解问题化为代数方程组，通过解代数方程组，得到出定解问题的解在离散点上的近似值组成的离散解，应用插值方法可从离散解得到定解问题在整个定解区域上的近似解。由此可见，有限差分方法有大体固定的模式，它有较强的通用性。但是，不能误认为不去了解这种逼近方法的基本知识，只是单纯模仿，便能轻易获得满意的结果。因为在应用这种逼近方法时会发生许多重要的但有时还是相当困难的数学问题，包括精度、稳定性与收敛性等。本章的目的是介绍有限差分方法的一些基本概念和构造差分格式的基本方法，以及给出一些力学问题有限差分求解的全过程。§1有限差分的基本概念在本节中，我们将首先针对单变量系统，基于Taylor展式导出若干有限差分表达式。接着，将所得结果推广到两个自变量以及多自变量函数系统。§1.1函数的表示首先考虑单变量函数)(xu，x为自变量。我们将区域x离散化为一系列点（或节点，有时也称作结点）（如图1.1），使得rrurhuxu==)()(，（,2,1,0=r）（1.1）图1.1等格距h时函数)(xu的离散化在用rh代替rx后，节点坐标就仅仅由整数r和格距h的乘积给定(这里不妨假设h为常数并规范化为小于1)，其中h称为沿x方向的步长。整数r表示节点沿x坐标的位置，通常，当0=x时0=r。当h是常数时，)(rhu就可表示为ru。(a)(b)图1.2（a）二维有限差分离散网格；（b）三维有限差分离散网格在二维情形下(图1.2a)，函数),(yxu在任何节点位置的值可表示表示为xu)(xuhh2h3h4h5rh)()(rhuxuurr==srsruskrhuyxu,),(),(==，),2,1,0;,2,1,0(==sr（1.2）其中x方向的格距为h，y方向的格距为k（或者，h称为沿x方向的步长，k称为沿y方向的步长）；整数r和s分别表示函数),(yxu沿x和y坐标的位置。此外，对于与任意点),(sr相邻的节点，我们可以表示如下],)1[(,1skhruusr±=±或])1(,[1,ksrhuusr±=±（1.3）类似地，同样可以对三维情形下（图1.2b)函数),,(zyxu给出离散表示，在此不再详述。§1.2单变量函数的有限差分公式Taylor级数展开式对于有限差分公式的建立和分类具有十分重要的意义，而且一定程度上其对于许多函数的近似逼近也是如此。因此，在具体研究各种数值方法之前，有必要对有关Taylor级数的知识作一下简单回顾或阐述以方便对此不甚熟悉的读者。单变量函数)(xu在离散点rx处的Taylor级数展开用我们前面所采用的记号可表示为++++=+rrrxxxxxxxxxrruhuhuhxuhxu!3!2)()(32（1.4a）或+−+−=−rrrxxxxxxxxxrruhuhuhxuhxu!3!2)()(32（1.4b）重新整理上两式后，可另写成−−−−+=rrrxxxxxxxrrxxuhuhhxuhxuu!3!2)()(2（1.5a）+−+−−=rrrxxxxxxxrrxxuhuhhhxuxuu!3!2)()(2（1.5b）进而，在点rx处的一阶导数的两个近似公式可由（1.5）给出huuhxuhxuuurrrrrxxxr−=−+≈=+1)()()(（1.6a）huuhhxuxuuurrrrrxxxr1)()()(−−=−−≈=（1.6b）由于上式中级数被截断，因此这些近似公式肯定存在一定的误差rE。此截断误差可由级数被截部分的第一项（也是昀大一项）表示出，即⎪⎪⎩⎪⎪⎨⎧−=+=−===)~(,)(2)~(,)(2~~rrxxxxrrxxxxrxxhxhouhhxxxhouhE（1.7）我们称此误差)(ho与h同阶。对于足够小的步长h，误差)(ho的绝对值将小于Ah（A为一常数）。（1.6a）和（1.6b）分别称作函数)(xu关于自变量x的一阶向前差商和向后差商。如果将（1.5a）与（1.5b）相加并求解rxu)(，则得huuurrrx2)(11−+−=（1.8）其被截去的第一项为)~(,)(6112~2+−==−=rrxxxxxrxxxhouhE（1.9）即截断误差为)(2ho阶。（1.8）称作函数)(xu关于自变量x的一阶中心差商。为了获得函数高阶导数的近似，例如二阶导数，我们将（1.5a）减去（1.5b）并求解rxxu)(，得到2112)(huuuurrrrxx−++−=（1.10）其被截去的第一项为)~(,)(12112~2+−==−=rrxxxxxxrxxxhouhE（1.11）即式（1.10）的截断误差为)(2ho阶。依次类推，我们不难得到更高阶（如：三阶、四阶）导数的近似表示公式。虽然我们可以沿用上面的方式推导出更为复杂的公式，但运算过程十分复杂，这里我们介绍另外一种使用算子推导的方法。首先，定义如下的符号以及算子（如表1.1）:表1.1差分算子以及符号符号算子差分表示式∆向前差分rrruuu−=∆+1∇向后差分1−−=∇rrruuuδ中心差分2/12/1−+−=rrruuuδE移位1+=rruEuµ平均2/)(2/12/1−++=rrruuuµD微分rxxxruxuDur)()d/d(===利用所定义的算子，可以用一种简单的推导和表述方式将所有可能的微分式的有限差分形式表示出来，其所得结果与通过Taylor级数所导出的公式相同。由上表中定义的各种线性算子，在不同的算子间存在很多关系式，这里给出几个作为例证。1)1(1−=∆⇒−=−=−=∆+EuEuEuuuurrrrrr（1.12）11111)1(−−−−−=∇⇒−=−=−=∇EuEuEuuuurrrrrr（1.13）22)(112/12/1−+−+−=+=rrrrruuuuuδδδµ（1.14）112/12/122)()(−+−++−=−==rrrrrrruuuuuuuδδδδ（1.15）22221123−−++−+−=rrrrruuuuuµδ（1.16）进一步，我们可将Taylor级数（1.4a）另表示为rrrxxxrxxrxrruhDuDhDhhDuhuhuhuu)exp()!3!21()(!3)(!2)(3322321=++++=++++=+（1.17）但由前面的定义，1+=rruEu，故有)exp(hDE=（1.18）或者⎪⎪⎩⎪⎪⎨⎧+∇+∇+∇=∇−−−∆+∆−∆=∆+==32323121)1ln(3121)1ln(lnEhD（1.19）同样地，容易证明)2sinh(2hD=δ（1.20）_____________________其证明如下：rrrxxxrxxrxrruhDuDhDhDhuhuhuhuu)2exp()!3)2/(!2)2/()2/(1()(!3)2/()(!2)2/())(2/(3322322/1=++++=++++=+（*1）rrrxxxrxxrxrruhDuDhDhDhuhuhuhuu)2exp()!3)2/(!2)2/()2/(1()(!3)2/()(!2)2/())(2/(3322322/1−=+−+−+−+=+−+−+−+=−（*2）由（*1）式与（*2）相减得到rrrrruuhDuhDuuδ=−−=−−+)2exp()2exp(2/12/1（*3）即有)2exp()2exp(hDhD−−=δ或)2sinh(2hD=δ（*4）获证。#___________________由式（1.20）,我们有−+−==−542321!523!321)2(sinh2δδδδhD（1.21）结合（1.19）式，可获得hD的各种不同表述，即函数一阶导数的几种差分表示式。在对这些公式进行平方、立方、四次方等，就可得到,,,443322DhDhDh。再将这些算子应用于ru或者2/1+ru，可得到单变量函数的前几阶导数表示式，由于高次方的运算相当冗长，这里我们只是给出如下的昀终结果。一阶导数：rrxuuh)3121()(32−∆+∆−∆=（1.22a）或rrxuuh)3121()(32+∇+∇+∇=（1.22b）或rrxuuh)301!31()(53−+−=δδδµ（1.22c）二阶导数：rrxxuuh)1211()(4322−∆+∆−∆=（1.23a）或rrxxuuh)1211()(4322+∇+∇+∇=（1.23b）或rrxxuuh)901121()(6422−+−=δδδ（1.23c）三阶导数：rrxxxuuh)4723()(5433−∆+∆−∆=（1.24a）或rrxxxuuh)4723()(5433+∇+∇+∇=（1.24b）或rrxxxuuh)120741()(7533−+−=δδδµ（1.24c）四阶导数：rrxxxxuuh)6172()(6544−∆+∆−∆=（1.25a）或rrxxxxuuh)6172()(6544+∇+∇+∇=（1.25b）rrxxxxuuh)240761()(8644−+−=δδδ（1.25c）从以上,)(,)(,)(32rxxxrxxrxuhuhuh的表达式中适当地截断无穷级数就可得到各种近似差分公式。例如，在向前差分公式（1.22a）中，若将一阶差分以后的各项舍去则得huuhuurrrrx−=∆=+1)(（1.26）上式与直接由Taylor级数方法得到的（1.6a）结果一致。其近似公式的截断误差由被略去的第一个差分项及其非零系数确定，即)~(,)()(21211~22+=−=∆−=rrxxxrrxxxhouhhuhE（1.27）这与Taylor级数误差估计结果也完全一致。表1.2中给出了单变量函数常用的前几阶有限差分近似表达公式以及相应的误差阶。表1.2单变量函数差分公式导数有限差分近似公式误差阶rxu)(huurr)(1−+)(hohuurr)(1−−)(hohuurr2)(11−+−)(2hohuuurrr2)34(12−+−++)(2hohuuuurrrr12)88(2112−−+++−+−4()ohrxxu)(211)2(huuurrr−++−)(2ho2211212)163016(huuuuurrrrr−−++−+−+−)(4horxxxu)(321122)22(huuuurrrr−−++−+−)(2horxxxxu)(42112)464(huuuuurrrrr−−+++−+−)(2ho§1.3多变量函数的有限差分公式利用上节的结果，我们可以直接推广到多变量函数),,,(321xxxu的许多有用的有限差分近似公式。这里，不妨以二元函数),(yxu为例，给出其多阶偏导数的差分近似表示公式。与单变量函数情形不同，需要补充的新概念包括了用符号sru,代替ru表示在节点),(sr处的函数值，以及对x求偏导数时y保持不变，反之亦然。结合huuurrrx)()(1−≈+以及图1.2a，不难得到)()(),(,,1,hohuuuxyxusrsrsrxsr+−==∂∂+（1.28a）)()(),(,1,,kokuuuyyxusrsrsrysr+−==∂∂+（1.28b）