-
抛物型方程是偏微分方程中的一种重要方程,它广泛应用于力学、天文学、物理学、生态学及工程技术等各个领域中.抛物型方程常见的数值解法有有限元法、有限体积法、边界元法和有限差分法等,在这些方法中,有限差分法仍然是求解抛物型方程的重要方法.有限差分法中常见的差分格式有古典显式格式、古典隐式格式和Crank-Nicolson格式等.古典显式格式的稳定性条件为r≤
$\frac{1}{2}$ ,截断误差仅为O(τ+h2).古典隐式格式和Crank-Nicolson格式则为无条件稳定的,古典隐式格式的截断误差为O(τ+h2),而Crank-Nicolson格式的截断误差为O(τ2+h2)[1-2].这些格式的精度都不高.因此,寻找稳定性好且精度高的差分格式就成了当前许多学者所研究的问题.本文主要考虑如下一维热传导方程
其中:u=u(x,t)表示温度,它是时间变量t与空间变量x的函数,a>0为热扩散项系数,它由材料的热传导率、密度与热容所决定,L是一个大于零的常量.为了确定一个具体的热传导过程,还必须知道物体的初始温度(初始条件),初始条件假设为
φ(x)是已知的光滑函数.另外物体在边界上所受到的外界的影响(边界条件)假设为
热传导方程除了用来刻画热量的传导过程,在实际问题中还有着广泛的应用.如生态学中提出的各种数学模型、神经轴突中电脉冲的传导及燃烧理论等许多物理现象,流体在多孔介质中的运动规律等都可归结为以上的热传导方程.对方程(1)的求解,有学者提出了改进的差分格式,精度比常见的经典格式精度要高[3-9],如:文献[8]给出了一族高精度恒稳格式,格式的截断误差达O(τ2+h6);文献[9]构造了一个截断误差达O(τ4+h4)的高精度隐格式,稳定性条件为0 < r≤0.5.本文构造了一个更高精度的三层九点隐式格式,格式的截断误差达到了O(τ3+h6),稳定性条件为r <
$\frac{\sqrt{5}}{10}$ .
全文HTML
-
设时间步长为τ,空间步长为h,网格区域由点集(xj,tn)(j=0,1,2,…,M;n=0,1,…)所组成,其中xj=jh,tn=nτ,h=
$\frac{L}{M}$ .并记ujn=u(xj,tn).用如下含参数的差分方程逼近微分方程(1)
其中:
${{\mathit{\Delta }}_{t}}u_{j}^{n}=\frac{u_{j}^{n+1}-u_{j}^{n}}{\tau }$ ,$\delta _{x}^{2}u_{j}^{n}=\frac{u_{j+1}^{n}-2u_{j}^{n}+u_{j-1}^{n}}{{{h}^{2}}}$ ,其余类推;ci(i=1,…,5)为待定系数.将(4)式中各节点上u的值在节点(jh,nτ)处作Taylor展开,并整理可得
为使格式(4)的截断误差达到O(τ3+h6),需下面方程组成立
在方程组(5)中,r=
$\frac{a\tau }{{{h}^{2}}}$ ,令c5=θ,可解得将所得各值代入(4)式,可得截断误差为O(τ3+h6)的一个隐式格式
-
利用Fourier分析法,可算出格式(6)的传播矩阵为
其中
传播矩阵G (s)的特征方程为
引理 1[10] 特征方程(7)的根满足|λ1,2|≤1的充要条件是
引理 2[10] 差分格式(6)稳定,即矩阵族Gn(s)(s∈[0, 1],n=1,2,…)一致有界的充要条件是
1) |λ1,2|≤1(λ1,2是方程(7)的两个根).
2) 使1-
$\frac{g_{11}^{2}}{4}$ =g112+4g12=0成立的s或不存在,或不属于区间[0, 1].定理 1 当0 < r <
$\frac{\sqrt{5}}{10}$ 且θ≥$\frac{40{{r}^{2}}-30r+1}{240{{r}^{2}}}$ 时,差分格式(6)稳定且收敛.证 首先考虑条件(2),当g12≠-1时,使1-
$\frac{g_{11}^{2}}{4}$ =g112+4g12=0的s不存在.再由条件(1)和式(8)知,格式(6)稳定的条件为-1+g12≤g11≤1-g12 < 2.由g11≤1-g12得
为确定起见,不妨假定分母
(10) 式成立的一个充分条件是
由(11),(12)式解得
当r <
$\frac{\sqrt{5}}{10}$ 时,$\frac{100{{r}^{2}}-60r+1}{240{{r}^{2}}}>\frac{480{{r}^{3}}-100{{r}^{2}}-108r-1}{480{{r}^{2}}\left( 1+3r \right)}$ ,故r <$\frac{\sqrt{5}}{10}$ 时(13)式优于(14)式,(13)式成立时(14)式也成立.而当(10)式成立时,由(9)式解得
(15) 式成立的条件为
又由1-g12 < 2可得
(16) 式成立的一个充分条件是
当r <
$\frac{1}{2}$ 时(17)式成立.由(18)式解得
而当r <
$\frac{1}{3}$ 时,$\frac{100{{r}^{2}}-60r+1}{240{{r}^{2}}}>\frac{40{{r}^{2}}-10r-9}{240{{r}^{2}}}$ 成立,故(13)式优于(19)式,(13)式成立时(19)式也成立.再由-1+g12≤g11得
(20) 式成立的一个充分条件是
由(21),(22)式解得
当r <
$\frac{1}{2}$ 时,$\frac{100{{r}^{2}}-60r+1}{240{{r}^{2}}}<\frac{40{{r}^{2}}-30r+1}{240{{r}^{2}}}$ 成立,故(23)式优于(13)式,(23)式成立时(13)式也成立.另外当r <
$\frac{\sqrt{5}}{10}$ 时,$\frac{140{{r}^{2}}-54r-120{{r}^{2}}-1}{480{{r}^{2}}}<\frac{40{{r}^{2}}-30r+1}{240{{r}^{2}}}$ ,故(23)式成立时(24)式也成立.综上所述,由Lax的稳定性与收敛性定理可得定理1结论成立.
特别地,当
$\theta =\frac{40{{r}^{2}}-30r+1}{240{{r}^{2}}}$ 时,差分格式(6)演化为
-
考虑扩散方程
利用格式(25)求数值解,并与精确解进行比较.
取h=
$\frac{1}{20}$ ,τ=rh2=$\frac{r}{400}$ ,r=$\frac{1}{10}$ ,$\frac{1}{8}$ ,$\frac{1}{6}$ ,$\frac{1}{5}$ 计算.为方便起见,用式(26)的精确解u(x,t)=e-tsinx计算第一层的值uj1.将本文格式(25)与文献[9]格式计算出的数值解和精确解进行比较,计算到n=200时的结果如表 1.由表 1看出,本文所构造差分格式的解与精确解有很好的吻合,与文献[9]格式相比,本文格式的精度比文献[9]格式的误差精度提高了10-3~10-4,结果说明本文构造的格式是一族高精度的差分格式.