求解三維空間分?jǐn)?shù)階對流擴(kuò)散方程的Douglas-Gunn格式

， ,，

(1.西北工業(yè)大學(xué) 計(jì)算科學(xué)研究中心 陜西 西安 710129； 2. 河南工業(yè)大學(xué) 理學(xué)院 河南 鄭州 450001)

0 引言

近年來，自然界中的反常擴(kuò)散現(xiàn)象受到科研人員的廣泛關(guān)注，為研究其獨(dú)特的物理過程，常常用分?jǐn)?shù)階偏微分方程建立相應(yīng)的數(shù)學(xué)模型. 其中，在包含對流和超擴(kuò)散兩個(gè)物理過程的散布現(xiàn)象中，粒子束的傳播與經(jīng)典的布朗運(yùn)動模型不再一致，此時(shí)把經(jīng)典對流擴(kuò)散方程中的空間二階導(dǎo)數(shù)替換成分?jǐn)?shù)階導(dǎo)數(shù)構(gòu)建的空間分?jǐn)?shù)階對流擴(kuò)散方程(SFADE)能更準(zhǔn)確地模擬這一現(xiàn)象. 分?jǐn)?shù)階導(dǎo)數(shù)或積分具有非局部性，這給相應(yīng)方程的求解帶來了很大困難，在大多數(shù)情況下很難得到解析解，因此研究可靠而有效的數(shù)值方法就顯得尤為重要. 目前常用的數(shù)值解法包括有限差分法[1-2]、有限元法[3-4]、有限體積法[5]、配點(diǎn)法[6]以及譜方法[7]等.

由于三維模型在科學(xué)研究中有廣泛應(yīng)用，本文考慮有限區(qū)域上帶有零Dirichlet邊界條件的三維SFADE的數(shù)值求解. 分?jǐn)?shù)階導(dǎo)數(shù)是一個(gè)非局部算子，這就使得離散SFADE得到的線性系統(tǒng)的剛度矩陣不再是稀疏矩陣，導(dǎo)致計(jì)算工作量和存儲量都非常大，尤其在多維空間情形下. 目前求解三維SFADE的數(shù)值方法還比較少. 文獻(xiàn)[8]采用了一種交替方向穩(wěn)定法(alternating direction implicit method, ADI)差分格式求解三維SFADE，并提高了精度.文獻(xiàn)[9]提出了一種求解三維分?jǐn)?shù)階擴(kuò)散方程的ADI差分格式.文獻(xiàn)[10]研究了一種求解三維空間分?jǐn)?shù)階擴(kuò)散方程的快速迭代ADI有限差分方法. 在本文中，我們將提出一種求解三維SFADE的有效的ADI有限差分格式，這種方法是將經(jīng)典的Douglas-Gunn格式中的二階中心差分算子推廣為包含左、右分?jǐn)?shù)階導(dǎo)數(shù)離散算子及一階中心差分算子在內(nèi)的復(fù)雜算子得到的，同時(shí)給出了該格式的穩(wěn)定性和收斂性的必要證明. 最后用數(shù)值實(shí)驗(yàn)驗(yàn)證了理論分析的結(jié)果.

1 三維SFADE及其Douglas-Gunn格式

本文考慮三維SFADE及它的初邊值條件為

(1)

u(x,y,z,t)=0, (x,y,z,t)∈?Ω×[0,T],

(2)

u(x,y,z,0)=u0(x,y,z), (x,y,z)∈Ω,

(3)

1.1 Riemann-Liouville分?jǐn)?shù)階導(dǎo)數(shù)的離散

其中系數(shù)為

記

(4)

(5)

同時(shí)，用中心差分近似對流項(xiàng)的一階導(dǎo)數(shù)，記

/2h1.

(6)

1.2 三維SFADE的有限差分近似

用公式(4)～(6)離散空間導(dǎo)數(shù)，時(shí)間方向采用Crank-Nicolson格式. 記

然后方程(1)就可以表示為

(7)

(8)

(9)

采用經(jīng)典的Douglas-Gunn格式分解式(9)得到ADI格式，即

(10)

(11)

(12)

(13)

與經(jīng)典的Douglas-Gunn格式不同，3個(gè)方向的二階中心差分算子在此處分別被替換為δα,x、δβ,y、δγ,z，它們是包含了左、右分?jǐn)?shù)階導(dǎo)數(shù)離散算子等在內(nèi)的復(fù)雜算子，可以認(rèn)為是經(jīng)典Douglas-Gunn格式在求解分?jǐn)?shù)階方程中的推廣. 接下來，我們將給出收斂性和穩(wěn)定性的必要證明.

2 收斂性和穩(wěn)定性分析

顯然，如果在格式(10)～(13)中消去中間解變量，則得到格式(9)，即格式(10)～(13)和(9)是等價(jià)的. 下面用矩陣法證明格式(9)是無條件穩(wěn)定和收斂的. 首先把方程(9)表示成矩陣形式， 令

(14)

(15)

(16)

(17)

其中Aα，Aβ，Aγ是Toeplitz矩陣，Aα和B分別表示為

其中：I是單位矩陣；符號?表示Kronecker積[12].Aβ、Aγ與Aα類似. 利用上述記號，式(9)可以寫為

τFn+1/2.

(18)

為了證明式(18)的穩(wěn)定性和收斂性，下面列出一些相關(guān)的引理和定理.

引理1[13]一個(gè)n階實(shí)矩陣A是正定的，當(dāng)且僅當(dāng)矩陣H=(A+AT)/2是正定的；H是正定的，當(dāng)且僅當(dāng)H的特征值都是正的.

引理2[13]設(shè)A是一個(gè)n階復(fù)矩陣，AH表示A的共軛轉(zhuǎn)置，記H=(A+AH)/2，則對于A的任意特征值λ，它的實(shí)部滿足不等式λmin(H)≤R(λ(A))≤λmax(H)，這里λmin(H)和λmax(H)分別表示H的最小和最大特征值.

煤炭資源影響著我國社會經(jīng)濟(jì)的發(fā)展，是社會發(fā)展的重要基礎(chǔ)資源。隨著我國工業(yè)技術(shù)的迅猛發(fā)展，煤炭資源的需求量逐漸增加，為新時(shí)代背景下煤炭資源的開發(fā)提出了更高的要求。隨著自動化技術(shù)在各領(lǐng)域中的廣泛應(yīng)用，煤礦機(jī)電設(shè)備的自動化程度也越來越高，有效滿足了煤炭資源供給需求。

引理3設(shè)A∈Rm×n,B∈Rr×s,C∈Rp×q，則(A?B)?C=A?(B?C).

證明此結(jié)論可以由Kronecker積的定義直接得到.

引理4設(shè)A,B∈Rm×n,C∈Rs×t，則有(A+B)?C=A?C+B?C，C?(A+B)=C?A+C?B.

證明此結(jié)論可以由Kronecker積的定義直接得到.

引理5[12]設(shè)A∈Rm×n,B∈Rr×s,C∈Rn×p，D∈Rs×t，則(A?B)(C?D)=AC?BD(∈Rmr×pt).

引理6[12]對于任意的矩陣A和B，有(A?B)T=AT?BT.

為了敘述并證明下述引理和定理，記‖·‖表示矩陣的2-范數(shù).

引理8設(shè)A∈Rn×n是正定矩陣，則對任意的σ∈R且σ>0，有‖(I+σA)-1‖<1.

引理9設(shè)A∈Rn×n是正定矩陣，則對任意的σ∈R且σ>0，有‖(I+σA)-1(I-σA)‖<1.

證明由矩陣2-范數(shù)的定義并記y=(I+σA)-1x，可得

引理10設(shè)A,B,I∈Rn×n，A和B乘積可交換，且(I-A)-1、(I-B)-1存在，則(I+A)與(I-B)-1，(I-A)-1與(I-B)-1也是乘積可交換的.

證明首先，由AB=BA，不難驗(yàn)證(I±A)(I-B)=(I-B)(I±A). 所以有

(I+A)(I-B)-1=(I-B)-1(I-B)(I+A)(I-B)-1=(I-B)-1(I+A)(I-B)(I-B)-1=

(I-B)-1(I+A)，(I-A)-1(I-B)-1=((I-B)(I-A))-1=((I-A)(I-B))-1=(I-B)-1(I-A)-1.

定理3差分格式(9)是無條件穩(wěn)定的.

(19)

(20)

3 數(shù)值結(jié)果

下面，我們通過兩個(gè)數(shù)值算例驗(yàn)證本文所提出的數(shù)值格式的穩(wěn)定性和收斂階，也就是說格式是有效的，并在時(shí)間和空間方向都具有較高的二階收斂精度. 設(shè)U和Uh分別表示問題(1)～(3)的解析解和采用格式(10)～(13)得到的數(shù)值解，用離散的L∞和L2范數(shù)計(jì)算全局截?cái)嗾`差，即

算例1在問題(1)～(3)中，取Ω=(0,1)×(0,1)×(0,1)，T=1；對流和擴(kuò)散系數(shù)分別為k1=k2=k3=-1，a1=a2=b1=b2=c1=c2=1；初值取u0(x,y,z)=x3(1-x)3y3(1-y)3z3(1-z)3. 已知的解析解為u(x,y,z,t)=e-tx3(1-x)3y3(1-y)3z3(1-z)3，由以上條件容易算出f(x,y,z,t).

對常系數(shù)算例1取優(yōu)化的步長比例M=N1=N2=N3進(jìn)行測試. 表1列出的數(shù)值結(jié)果表明，用格式(10)～(13)計(jì)算常系數(shù)問題(1)～(3)時(shí)，算法是無條件穩(wěn)定的，而且在時(shí)間及空間方向都是二階收斂的，這和理論分析的結(jié)果一致.

算例2在問題(1)～(3)中，取Ω=(0,1)×(0,1)×(0,1)，T=1；對流和擴(kuò)散系數(shù)分別為k1=0.25x，k2=0.25y，k3=0.25z，a1=xα，a2=(1-x)α，b1=yβ，b2=(1-y)β，c1=zγ，c2=(1-z)γ；初值取u0(x,y,z)=x3(1-x)3y3(1-y)3z3(1-z)3．已知的解析解為u(x,y,z,t)=e-tx3(1-x)3y3(1-y)3z3·(1-z)3，由以上條件f(x,y,z,t)容易算出.

表2列出了變系數(shù)算例2的數(shù)值結(jié)果，這里也取優(yōu)化的步長比例M=N1=N2=N3進(jìn)行測試，數(shù)值結(jié)果表明用格式(10)～(13)計(jì)算變系數(shù)問題(1)～(3)時(shí)，算法是無條件穩(wěn)定的，而且在時(shí)間及空間方向也都具有二階收斂率，這和理論分析的結(jié)果是非常吻合的.

表1 算例1在時(shí)刻t=1，取M=N1=N2=N3的數(shù)值誤差和收斂階Tab.1 The errors and convergence orders for example 1 at t=1 with M=N1=N2=N3

表2 算例2在時(shí)刻t=1，取M=N1=N2=N3的數(shù)值誤差和收斂階Tab.2 The errors and convergence orders for example 2 at t=1 with M=N1=N2=N3

4 結(jié)論

本文將求解三維整數(shù)階拋物方程的經(jīng)典Douglas-Gunn格式推廣到分?jǐn)?shù)階，提出了一種求解三維SFADE的有效的數(shù)值方法，并證明該格式具有無條件穩(wěn)定性和較高的二階收斂精度，必要而充足的數(shù)值實(shí)驗(yàn)驗(yàn)證了理論結(jié)果. 最后，由于分?jǐn)?shù)階導(dǎo)數(shù)是非局部算子，對于多維空間問題的求解需要耗費(fèi)較大的計(jì)算工作量和空間存儲量，在今后的工作中，我們將考慮開展適當(dāng)?shù)目焖偎惴?，以減少計(jì)算花費(fèi)和加快計(jì)算速度.