有限差分方法
一种求偏微分(或常微分)方程和方程组定解问题的数值解的方法,简称差分方法。
微分方程的定解问题就是在满足某些定解条件下求微分方程的解。在空间区域的边界上要满足的定解条件称为边值条件。如果问题与时间有关,在初始时刻所要满足的定解条件,称为初值条件。不含时间而只带边值条件的定解问题,称为边值问题。与时间有关而只带初值条件的定解问题,称为初值问题。同时带有两种定解条件的问题,称为初值边值混合问题。
定解问题往往不具有解析解,或者其解析解不易计算。所以要采用可行的数值解法。有限差分方法就是一种数值解法,它的基本思想是先把问题的定义域进行网格剖分,然后在网格点上,按适当的数值微分公式把定解问题中的微商换成差商,从而把原问题离散化为差分格式,进而求出数值解。此外,还要研究差分格式的解的存在性和唯一性、解的求法、解法的数值稳定性、差分格式的解与原定解问题的真解的误差估计、差分格式的解当网格大小趋于零时是否趋于真解(即收敛性),等等。
有限差分方法具有简单、灵活以及通用性强等特点,容易在计算机上实现。
偏微分方程初值问题的差分法 许多物理现象随着时间而发生变化、如热传导过程、气体扩散过程和波的传播过程都与时间有关。描述这些过程的偏微分方程具有这样的性质:若初始时刻t=t的解已给定,则t>t时刻的解完全取决于初始条件和某些边界条件。利用差分法解这类问题,就是从初始值出发,通过差分格式沿时间增加的方向,逐步求出微分方程的近似解。
双曲型方程的差分方法 最简单的双曲型方程的初值问题是:
式中(x)为已知初值函数。这初值问题的解是:
u(x,t)=(x-at)。(2)
由(2)可见,(1a)(1b)的解(2)当a>0时代表一个以有限的速度a沿特征线x-at=常数向右传播的波,而解u(x,t)在P(,)点的值完全由(x)在x轴上的点A(-а,0)的值决定。A点就是双曲型方程(1a)在P点的依赖域(图1)。现以初值问题(1)为例介绍初值问题差分方法的基本思想。①剖分网格 用网格覆盖(1a),(1b)的定解区域,如图2所示,在x,t平面的上半部作两族平行于坐标轴的直线:
x=x=jΔx,
j=0,±1,±2,…
t=t=nΔt,
n=0,1,2,…
并称之为网格线。Δx,Δt分别称为空间步长和时间步长。网格线的交点(jΔx,nΔt)称为格点。②建立差分格式 以下除特别声明外,总设a>0,由泰勒公式,有:
,(3a)
和
-1<θ<0(3b)
解出,代入(1a)得:
(4)
式中
(5)
E是微分方程(1a)用它的解在相邻三个格点(见图2)上的值的差分来表示的形式。略去(4)中关于Δx,Δt的高阶项E,得到一个较简单的差分方程,但微分方程的解u(jΔx,nΔt)不再是这方程的解;设这个方程的解是u,u满足的方程是:
;(6)
式(6)还可写成:
;(6′)
初值条件(1b)此时就是:
u=(x),j=0,±1,±2,…。(7)
差分方程(6)和相应的初值条件(7)合称差分格式,利用这些格式可逐步算出t=Δt,2Δt,…各时间层的u,u,…,等等。这个把微分方程化为近似的差分方程的过程常称为离散化。
③差分格式的截断误差和相容性 (5)中的E是把微分方程充分光滑的解代入差分方程(6)的结果,它说明微分方程(1a)和差分方程(6)的区别,称为差分格式(6)的截断误差。式(6)的截断误差对Δt和Δx都是一阶的,写成O(Δx+Δt),因此称差分格式(6)为一阶相容格式。一般说,如果Δx,Δt趋于零,截断误差也趋于零,则差分方程与微分方程是相容的。不相容的格式的解不能作为原微分方程的近似解,因而是无用的。方程(1a)的离散化过程也不是唯一的。例如取数值微分公式:
代替微分方程(1a)中的,可得另一个差分方程:
,(9)
它的截断误差是O(Δx+Δt)阶的,也是相容的差分格式,再若用数值微分公式
代替(1a)中的,又得到截断误差为O(Δx+Δt)的相容差分格式:
。(11)
但是,并不是每个相容格式都有用。
④差分格式的收敛性 设P(,)是求解区域中的一点,取步长Δx,Δt使=jΔx,=nΔt,用差分格式算出u,如果当Δx,Δt→0时u-u(,)→0,便可用步长充分小时的u作为微分方程的解u(jΔx,nΔt)的近似,这种差分格式便是收敛的。
双曲型微分方程的解,对求解区域内一点(,)而言,在初值区域内有一个依赖域,差分方程也是如此,对于差分方程(6),点(jΔx,nΔt)的依赖域是初值线上区间[(j-n)Δx,jΔx]。如令Δt/Δx=r=常数,=jΔx,=nΔt,则差分方程(6)在点(,)的依赖域为[-a/r,],并且步长比r固定时,依赖域与Δx,Δt无关。
差分方程(9)在(,)的依赖域是[-a/r,+a/r],而差分方程(11)的依赖域则是[,+a/r],R.库朗等人曾经证明,差分格式收敛的一个必要条件是差分方程的依赖域应包含微分方程的依赖域,这个条件叫作“库朗条件”。从图3中可以看到,对于差分方程(6),这个条件是-a/r≤-a≤,即对于格式(9),库朗条件是,两者不同。对于格式(11),库朗条件是≤-a≤+a/r;在a>0时,显然不能成立,所以格式(11)当a>0时不收敛,因而也是无用的。格式(6)a>0在而库朗条件满足时,的确是收敛的。因为的离散化误差
适合由此可知:
又因差分格式与微分方程的初值相同,于是可知:
这说明条件
满足时,格式(6)收敛。如果a<0,格式(6)不收敛。但当时,
格式(11)收敛。这两个格式称为“迎风格式”,因为a>0时,用向后差商代替,往上风取近似值;当a<0时则用向前差商代替,也是往上风取近似值。可见作(1)的差分格式时,要考虑波的传播方向。
⑤差分格式的稳定性 用一个差分格式计算时,初值的误差必然要影响到以后各层。通常希望这误差的影响不会越来越大,以致完全歪曲了差分方法的真解,这便是稳定性问题。讨论时,常把问题化简,设初值有误差ε,而以后的计算并不产生误差,由于误差ε,使变成了+ε,但+ε仍满足所适合的差分格式。定义一种衡量 t=t层格点上ε的大小的所谓范数‖ε‖,若有常数K>0 使当 Δt、Δx→0 而0≤t=nΔt≤T 时,恒有‖ε‖≤K‖ε‖,则称此差分格式是稳定的。以格式(6)为例,ε适合差分方程:
可证当时,取,则有
这说明,用格式(6)计算时,若步长比合于库朗条件,则初值误差的影响不增长,即使Δt缩小,算到t=T时,也不再增大,因而格式是稳定的。
对于线性偏微分方程组的稳定性理论,J.von诺伊曼曾用傅里叶分析作了系统研究,把差分方程的解表成谐波的叠加,考察其中一个谐波
(12)
的增长情况,式中k为实数;G=G(k,Δt)称为增长因子。若对于一切谐波,(12)的振幅一致有界,即对一切合于0≤nΔt≤T的n和充分小的Δt都有|G|≤K,K为常数,则此差分格式是稳定的。具体地说,对格式(6),把(12)代入(6),得:
而
故当时,G≤1,解的振幅不增加,所以格式(6)是稳定的。
相容性和库朗条件都不能保证稳定性,例如对格式(9),把(12)代入,得:
而
故当sinkΔx0时,恒有|G|>1,解的振幅逐层增加,所以虽然格式(9)是相容的格式,并且适合库朗条件,但它仍是不稳定的,因而也是无用的。
P.D.拉克斯1956年曾证明:对于线性偏微分方程组的适定的初值问题,一个与之相容的线性差分格式是收敛的格式的充分必要条件是这格式的稳定性。
非线性问题没有相应的等价定理。
抛物型方程的差分方法 抛物型方程的定解问题是初值问题或初值边值问题。为了说明抛物型方程差分方法的某些特点,考虑热传导方程的初值、边值问题:
式中a>0为常数;(x)为给定的连续函数。这里也是用直线
x=x=jΔx,j=0,1,…,Μ,
t=t=nΔt,n=0,1,2,…,N(=T/Δt),
剖分求解区域为矩形网格(见图4),式中Δx=1/Μ,Μ为正整数。利用数值微分公式:
中的(14a)及(14b),从微分方程(13)可得差分格式:
这里计算u时只用到前一层的u,u及u,是一个显示格式。要是用(14a)和(14c),则得到一个隐式差分格式:
此时计算u必须解一个线性代数方程组。用诺伊曼方法,可以证明:当常数时,对于任何r值,格式(16)是无条件稳定的;格式(15)则是条件稳定的,稳定性条件是。显然格式(15)的时间步长要受 到较大的限制。隐式格式归结为解线性代数方程组。最简便的方法是追赶法,其要点如下:对于每个时间步n+1,把(16)改写成:
用消去法可导出如下两套递推公式:
和
用公式(18a)从j=0出发,逐步先求出α和β,然后用(18b)从j=Μ-1逐步求出u,u,…,u,u来。当A>0,B>0,C≥A+B时,这两个递推过程都是稳定的。
偏微分方程边值问题的差分法 物理上的定常问题,如弹性力学中的平衡问题,亚声速流、 不可压粘性流、电磁场及引力场等可归结为椭圆型方程。其定解问题为各种边值问题,即要求解在某个区域D内满足微分方程,在边界上满足给定的边界条件。椭圆型方程的差分解法可归结为选取合理的差分网格,建立差分格式,求解代数方程组以及考察差分格式的收敛性等问题。
泊松方程是椭圆型方程的典型例子,它的第一边值问题为:
式中D为x,y平面上某个封闭区域;дD为它的边界(图5);(x,y),g(x,y)为连续函数;u(x,y)为未知解。以下简单介绍其差分解法的基本思想。
通常可将定解区域剖分成矩形网格或三角形网格。三角形网格对不规则区域较为方便。为简便起见,设D为单位正方形,x 和y方向均取为等距步长h,并用直线x=ih,y=jh(i,j=0,1,2,…,N)将此正方形D={0≤x≤y≤1}剖分成正方形网格。
在格点(i,j)上,微商u、u分别用x、y方向的二阶中心差商来代替,得到差分格式:
的截断误差是二阶的。将u按顺序u,u,…,u,u,u,…,排列,并用矢量形式表示为:
u=(u,u,…,u,u,u,…,u)T代表转置,则(20)可写成方程组
,(21)
其中系数矩阵A为三对角块状方阵
,I为(N-1)×(N-1)的单位矩阵,S亦为(N-1)×(N-1)矩阵。
偏微分方程边值问题的差分方程组的特点是系数矩阵中非零元素很少,即是稀疏矩阵。近年来由于稀疏矩阵技术的发展,解差分方程组时,直接法受到了较多的重视。迭代法是用逐次逼近的方式得到差分方程组的解,它的存储量小,程序简单,因此常用于椭圆型差分方程组的求解。迭代方法很多,最基本的有三种:
①同时位移法(也称雅可比法):
,(23)
n代表迭代的次数。
②逐个位移法(也称赛德耳法):
,(24)
当已算出时,取,否则。
③松弛法:
,(25)
当已算出时,,否则,
,
式中ω为松弛因子;0<ω<2,ω<1为低松弛法,ω>1为超松弛法。三个方法中超松弛法收敛最快,是常用的方法之一。
差分方法的发展和应用 前面阐述了两个自变量,线性方程的差分法。实际问题常会遇到多个自变量,非线性的方程或方程组;它们还可能是混合型的偏微分方程(如机翼的跨声速绕流),其解包含着各种间断(如激波间断、按触间断等)。非线性问题的差分法求解是十分困难的。随着电子计算机的发展,在解决各种非线性问题中,差分法得到了很快的发展,并且出现了许多新的思想和方法,如守恒差分格式,时间相关法、分步法等。
守恒差分格式 数学物理偏微分方程通常代表某种物理、力学中的守恒律(如质量守恒、动量守恒、能量守恒、粒子守恒等)。原始问题的差分格式,若能保持同样的守恒性质,就称为守恒差分格式。守恒性反映出物理问题的整体性质,用它来检验差分格式的好坏是合理的。对于间断的问题,守恒格式特别重要。从积分守恒关系式出发,利用积分插值方法容易得到守恒格式。这时对于复杂的求解区域、各种类型边界条件、间断系数等复杂情况都可以处理。
时间相关法 把定常的微分问题用一个相应的非定常问题来代替,然后用差分法解后者的初值问题,要求当t→∞时,它的稳定解为原来问题的解,这类方法叫作时间相关法。实践上,当计算时间足够大时,就能得到满足给定精度的近似解。例如拉普拉斯方程第一边值问题:
可以用热传导方程的初边值问题:
(27)
来代替。若用显式格式计算(27),可避免解大型代数方程组。特别是当微分方程的类型在定解区域内发生变化时,可只用一种类型来算,而使问题大大化简。这种方法在定常问题中广泛使用。缺点是达到定常解的计算时间较长,有待改进。
分步法 把复杂的问题的每一时间步分解成几个中间步,例如把多维问题按坐标分解为几个一维问题,然后用差分法解这些比较简单的各中间步,最后得到原始问题的近似解,这类方法叫作分步法。交替方向法、预估-修正法、时间分裂法、因式分解法等都属此类。以二维抛物型方程定解问题:
为例,用显式格式求解,时间步长受稳定性条件:
的限制,用隐式格式,则归结为大型线性代数方程组,解起来比较麻烦。1955年皮斯曼-拉什福德提出交替方向隐式格式:
(i=1,2,…,N-1;j=1,2,…,Μ-1;n=0,1,2,…)
δ为中心差分算符,第一步x方向取隐式,y方向取显式,第二步则相反。两步合成无条件稳定的格式。由于每一步可用追赶法求解,大大简化了解法。交替方向法出现后,进一步发展了各种形式的分步格式,并可推广到任何维数的方程或方程组的情形,困难在于边界条件的处理。
有限差分方法已成为解各类数学物理问题的主要数值方法,也是计算力学中的主要数值方法之一。有些解偏微分问题的方法(如特征线法、直线法)实质上也是差分方法的一种形式。在固体力学中,有限元方法出现以前,主要采取差分方法;在流体力学中,差分方法仍然是主要的数值方法。当然,对于某些具有复杂的几何形状及复杂的流动现象的实际问题,差分方法还有待进一步发展。
参考书目 冯康等编:《数值计算方法》,国防工业出版社,北京,1978。 胡祖炽编:《计算方法》,高等教育出版社,北京,1959。 清华大学、北京大学《计算方法》编写组编:《计算方法》,科学出版社,北京,1980。 朱幼兰等著:《初边值问题差分法及绕流》,科学出版社,北京,1980。 R.D.里奇特迈尔著,何旭初等译:《初值问题差分方法》,科学出版社,北京,1966。(R.D.Richtmyer,Difference Methods for Initial-Value Problems,Interscience Pub.,New York,1957.) R.D.Richtmyer,K.W.Morton, Difference Methods for Initial-Value Problems,2nd ed.,Interscience Pub.,New York,1967.