x=ih1,i=0,1,…,N1,h1=a/N1
y=jh2,j=0,1,…,N2,h2=b/N2
将G剖分为网格区域,见图9-1。
h1,h2分别称为x方向和y方向的剖分步长,网格交点(xi,yi)称为剖分节点(区域内节点集合记为Gh={(xi,yi);(xi,yi)∈G}),网格线与边界Γ的交点称为边界点,边界点集合记为Γh。
现在将微分方程(9.1)在每一个内节点(xi,yi)上进行离散。
在节点(xi,yi)处,方程(9.1)为
(9.5)
需进一步离散(9.5)中的二阶偏导数。
为简化记号,简记节点(xi,yi)=(i,j),节点函数值u(xi,yi)=u(i,j)。
利用一元函数的Taylor展开公式,推得二阶偏导数的差商表达式
代入(9.5)式中,得到方程(9.1)在节点(i,j)处的离散形式
其中。
舍去高阶小项,就导出了u(i,j)的近似值ui,j所满足的差分方程
(9.6)
在节点(i,j)处方程(9.6)逼近偏微分方程(9.1)的误差为,它关于剖分步长是二阶的。
这个误差称为差分方程逼近偏微分方程的截断误差,它的大小将影响近似解的精度。
在差分方程(9.6)中,每一个节点(i,j)处的方程仅涉及五个节点未知量ui,j,ui+1,j,ui-1,j,ui,j+1,ui,j-1,因此通常称(9.6)式为五点差分格式,当h1=h2=h时,它简化为
差分方程(9.6)中,方程个数等于内节点总数,但未知量除内节点值ui,j,(i,j)∈Gh外,还包括边界点值。
例如,点(1,j)处方程就含有边界点未知量u0,j。
因此,还要利用给定的边值条件补充上边界点未知量的方程。
对于第一边值条件式(9.2),可直接取ui,j=α(xi,yi),(i,j)∈Γh(9.7)
对于第三(k=0时为第二)边值条件式(9.4),
以左边界点(1,j)为例,见图9-2,
利用一阶差商公式
则得到边界点(0,j)处的差分方程
(9.8)
联立差分方程(9.6)与(9.7)或(9.8)就形成了求解Poisson方程边值问题的差分方程组,它实质上是一个关于未知量{ui,j}的线性代数方程组,可采用第2,3章介绍的方法进行求解。
这个方程组的解就称为偏微分方程的差分近似解,简称差分解。
考虑更一般形式的二阶椭圆型方程
(9.9)
其中A(x,y)≥Amin>0,B(x,y)≥Bmin>0,E(x,y)≥0。
引进半节点
利用一阶中心差商公式,在节点(i,j)处可有
对类似处理,就可推得求解方程(9.9)的差分方程
(9.10)
其中
(9.11)
显然,当系数函数A(x,y)=B(x,y)=1,C(x,y)=D(x,y)=E(x,y)=0时,椭圆型方程(9.9)就成为Poisson方程(9.1),而差分方程(9.10)就成为差分方程(9.6)。
容易看出,差分方程(9.10)的截断误差为阶。
9.1.2一般区域的边界条件处理
前面已假设G为矩形区域,现在考虑G为一般区域情形,这里主要涉及边界条件的处理。
考虑Poisson方程第一边值问题
(9.12)
其中G可为平面上一般区域,例如为曲边区域。
仍然用两组平行直线:
x=x0+ih1,y=y0+jh2,i,j=0,±1,…,对区域G进行矩形网格剖分,见图9-3。
如果一个内节点(i,j)的四个相邻节点(i+1,j),(i-1,j),(i,j+1)和(i,j-1)属于,则称其为正则内点,见图9-3中打“。
”号者;如果一个节点(i,j)属于且不为正则内点,则称其为非正则内点,见图9-3中打“.”号者。
记正则内点集合为,非正则内点集合为。
显然,当G为矩形区域时,成立。
在正则内点(i,j)处,完全同矩形区域情形,可建立五点差分格式
(9.13)
在方程(9.13)中,当(i,j)点临近边界时,将出现非正则内点上的未知量,因此必须补充非正则内点处的方程。
若非正则内点恰好是边界点,如图9-4中
D点,则利用边界条件可取uD=α(D)
对于不是边界点的非正则内点,如图9-4中B点,一般可采用如下两种处理方法。
a.直接转移法.取与点B距离最近的边界点(如图9-4中E点)上的u的值作为u(B)的近似值uB,即uB=u(E)=α(E)
直接转移法的优点是简单易行,但精度较低,只为一阶近似。
b.线性插值法.取B点的两个相邻点(如图9-4中边界点A和正则内点C作为插值节点对u(B)进行线性插值
则得到点B处的方程
线性插值法精度较高,为二阶近似。
对每一个非正则内点进行上述处理,将所得到的方程与(9.13)式联立,就组成了方程个数与未知量个数相一致的线性代数方程组。
求解此方程组就可得到一般区域上边值问题(9.12)的差分近似解。
对于一般区域上二阶椭圆型方程(9.9)的第一边值问题,可完全类似处理。
第二、三边值条件的处理较为复杂,这里不再讨论。
9.2抛物型方程的差分方法
本节介绍抛物型方程的差分方法,重点讨论差分格式的构造和稳定性分析。
9.2.1一维问题
作为模型,考虑一维热传导的初边值问题
(9.14)
(9.15)
(9.16)
其中a是正常数,都是已知的连续的函数。
现在讨论求解问题(9.14)-(9.18)的差分方法。
首先对求解区域G={0≤x≤l,0≤t≤T}进行网格剖分。
取空间步长h=l/N,时间步长τ=T/M,其中N,M是正整数,作两族平行直线
将区域G剖分成矩形网格,见图9-5,网格交点(xj,tk)称为节点。
用差分方法求解初边值问题(9.14)-(9.16)就是要求出精确解u(x,t)在每个节点(xj,tk)处的近似值。
为简化记号,简记节点(xj,tk)=u(j,k)。
利用一元函数的Taylor展开公式,可推出下列差商表达式
(9.17)
(9.18)
(9.19)
(9.20)
1.古典显格式
在区域G的内节点(j,k)处,利用公式(9.17)和(9.20),可将偏微分方程(9.14)离散为
其中。
舍去高阶小项,就得到节点近似值(差分解)所满足的差分方程
(9.21)
显然,在节点(j,k)处,差分方程(9.21)逼近偏微分方程(9.14)的误差为,这个误差称为截断误差,它反映了差分方程逼近偏微分方程的精度。
现将(9.21)式改写为便于计算的形式,并利用初边值条件(9.15)与(9.16)补充上初始值和边界点方程,则得到
(9.22)
其中称为网比。
与时间相关问题差分方程的求解通常是按时间方向逐层进行的。
对于差分方程(9.22),当第k层节点值已知时,可直接计算出第k+1层节点值。
这样,从第0层已知值开始,就可逐层求出各时间层的节点值。
差分方程(9.22)的求解计算是显式的,无须求解方程组,故称为古典显格式。
此外,在式(9.22)中,每个内节点处方程仅涉及k和k+1两层节点值,称这样的差分格式为双层格式。
差分方程(9.22)可表示为矩阵形式
(9.23)
其中
2.古典隐格式
在区域G的内节点(j,k)处,利用公式(9.18)和(9.20),可将偏微分方程(9.14)离散为
舍去高阶小项,则得到如下差分方程
(9.24)
它的截断误差为,逼近精度与古典显格式相同。
改写(9.24)式为便于计算的形式,并补充上初始值与边界点方程,则得到
(9.25)
与古典显格式不同,在差分方程(9.25)的求解中,当第k-1层值已知时,必须通过求解一个线性方程组才能求出第k层值,所以称(9.25)式为古典隐格式,它也是双层格式。
差分方程(9.25)的矩阵形式为
(9.26)
其中
向量同(9.23)式中定义。
从(9.26)式看到,古典隐格式在每一层计算时,都需求解一个三对角形线性方程组,可采用追赶法求解。
3.Crank-Nicolson格式(六点对称格式)
利用一元函数Taylor展开公式可得到如下等式
使用这两个公式,在点离散偏微分方程(9.14),然后利用(9.20)式进一步离散二阶偏导数,则可导出差分方程
(9.27)
其截断误差为,在时间方向的逼近阶较显格式和隐格式高出一阶。
这个差分格式称为Crank-Nicolson格式,有时也称为六点对称格式,它显然是双层隐式格式。
改写(9.27)式,并补充初始值和边界点方程得到
(9.28)
它的矩阵形式为
(9.29)
在每层计算时,仍需求解一个三对角形方程组。
4.Richardson格式
利用公式(9.19)和(9.20),可导出另一个截断误差为阶的差分方程
称之为Richardson格式。
可改写为
(9.32)
这是一个三层显式差分格式。
在逐层计算时,需用到和两层值才能得到k+1层值。
这样,从第0层已知值开始,还须补充上第一层值,才能逐层计算下去。
可采用前述的双层格式计算。
除上述四种差分格式外,还可构造出许多逼近偏微分方程(9.14)的差分格式,但并不是每个差分格式都是可用的。
一个有实用价值的差分格式应具有如下性质:
(1)收敛性。
对任意固定的节点(xj,tk),当剖分步长时,差分解应收敛到精确解u(xj,tk)。
(2)稳定性。
当某一时间层计算产生误差时,在以后各层的计算中,这些误差的传播积累是可控制而不是无限增长的。
理论上可以证明,在一定条件下,稳定的差分格式必然是收敛的。
因此,这里主要研究差分格式的稳定性。
作为例子,先考查Richardson格式的稳定性。
设是当计算过程中带有误差时,按Richardson格式(9.30)得到的实际计算值。
是理论值,误差。
假定右端项的计算是精确的,网比,则满足
(9.31)
设前k-1层计算时精确的,误差只是在第k层点发生,即
。
则利用(9.31)式可得到误差的传播情况,见表9-1。
表9-1r=1/2时Richardson格式的误差传播
j
k
j0-4
j0-3
j0-2
j0-1
j0
j0+1
j0+