陳景華，陳雪娟

(1.集美大學(xué)理學(xué)院，福建 廈門 361021；2.集美大學(xué)理學(xué)院數(shù)字福建大數(shù)據(jù)建模與智能計算研究所，福建 廈門 361021)

0 引言

分?jǐn)?shù)階微分方程是廣義的非整數(shù)階的微分方程。由于分?jǐn)?shù)階算子的非局部性，分?jǐn)?shù)階微分方程非常適合用來描述現(xiàn)實世界中具有記憶及遺傳性質(zhì)的材料[1-3]。 然而，還有很多問題不能采用分?jǐn)?shù)階系統(tǒng)來解決，如非均質(zhì)多孔和裂隙介質(zhì)中的溶質(zhì)遷移試驗[4-5]。由于復(fù)雜流體流動的力學(xué)性質(zhì)會隨時空尺度而發(fā)生一定的變化，因此污染物等溶質(zhì)擴散將產(chǎn)生所謂的多尺度效應(yīng)。 近年來，研究者們[6-9]采用Riesz分布階的分?jǐn)?shù)階擴散方程(其中分布階導(dǎo)數(shù)在數(shù)學(xué)上可以看成是對分?jǐn)?shù)階導(dǎo)數(shù)在給定范圍內(nèi)關(guān)于其階數(shù)的一個積分，即分?jǐn)?shù)階導(dǎo)數(shù)是在一個單位區(qū)間上的積分)來描述這種反常擴散。Diethelm等[10]研究了Caputo分布階常微分方程的數(shù)值解并進(jìn)行了收斂性分析；Atanackovic等[11]證明了時間分布階柯西問題解存在性，并求出解析解；Luchko等[12]考慮了一個開區(qū)間上的分布階時間分?jǐn)?shù)階擴散方程，并證明了這個問題的解的唯一性和存在性；Jiao等[9]給出了分布階導(dǎo)數(shù)在信號處理中的分布階的濾波器和控制系統(tǒng)中的最優(yōu)分布階阻尼兩個應(yīng)用。 Podlubny等[13]用矩陣方法離散分布階導(dǎo)數(shù)和積分,給出數(shù)值例子驗證算法的精確性，但是沒有具體的數(shù)值理論分析；Gorenflo等[14]研究了一種分布階時間分?jǐn)?shù)階擴散波動方程，并利用傅里葉變換和拉普拉斯變換技術(shù)得到了此問題的基本解?，F(xiàn)有文獻(xiàn)關(guān)于分布階微分方程數(shù)值方法的研究還比較少。2014年，Katsikadelis[15]研究了Caputo分布階線性與非線性常微分方程的數(shù)值解，但沒有進(jìn)行理論分析。2018年，Semarya等[16]運用移位的分?jǐn)?shù)階積分切比雪夫算子矩陣與修正的Picard數(shù)值方法求解分布階線性分?jǐn)?shù)階常微分方程，但沒有進(jìn)行理論分析。分布階微分算子作為分?jǐn)?shù)階微分算子的發(fā)展形式，已經(jīng)越來越受到研究人員的重視。相對于固定階導(dǎo)數(shù)，這類算子的應(yīng)用領(lǐng)域更為廣闊，但由于算子定義的復(fù)雜性，相應(yīng)的數(shù)值算法更為復(fù)雜，并不能簡單照搬分?jǐn)?shù)階的數(shù)值算法。本文考慮用有效的數(shù)值方法求解空間分布階的分?jǐn)?shù)階擴散方程。

考慮以下Riesz空間分布階的分?jǐn)?shù)階擴散方程[17-18]：

(1)

邊界條件：

u(xL,t)=0,u(xR,t)=0,t∈(0,T],

(2)

初始條件：

u(x,0)=φ(x),x∈[xL,xR],

(3)

這里P(α)是非負(fù)加權(quán)函數(shù)[10]，滿足：

(4)

(5)

這里cα=-1/(2cos(πα/2)),且

(7)

1 積分項的離散

首先，把積分區(qū)間劃分成若干個子區(qū)間：1=α0<α1<…<αS=2,S∈N，并記

Δαi=αi-αi-1=1/S=σ,αi-1/2=(αi-1+αi)/2,i=1,2,…,S。

(8)

(9)

(10)

則Riesz空間分布階的分?jǐn)?shù)階擴散方程可以離散成具有多項分?jǐn)?shù)階導(dǎo)數(shù)的微分方程：

?u(x,t)/?t=Dxu+f(x,t)+O(σ4)

(11)

具有邊界和初始條件(2)～(3)。

2 數(shù)值離散格式

首先對空間和時間進(jìn)行離散。取正整數(shù)N,設(shè)時間步長τ=T/N,tn=nτ,0≤n≤N。記時間區(qū)域Ωτ={tn|0≤n≤N}。取正整數(shù)M，并設(shè)空間步長h=(xR-xL)/M,xi=xL+ih,0≤i≤M?？臻g區(qū)域Ωh={xi|0≤i≤M}?？臻g平均差分算子Aα定義為：

平均算子A定義為：

(12)

(13)

(14)

(15)

(16)

其中,

(17)

(18)

且

(19)

(20)

及差分算子δxu：

(21)

(22)

其中,

(23)

由方程(9)～(16)得到：

A(Dxu(xi,t))=δxu(xi,t)+O(h4)。

(24)

考慮方程(1)在點(xi,t)處的方程為：

ut(xi,t)=Dxu(xi,t)+f(xi,t)+O(σ4),

(25)

在方程(25)兩端取平均算子A，則有：

Aut(xi,t)=ADxu(xi,t)+Af(xi,t)+O(σ4),xi∈Ωh。

(26)

由式(24)可得：

Aut(xi,t)=δxu(xi,t)+Af(xi,t)+O(σ4+h4),xi∈Ωh。

(27)

下面對時間離散。對于方程(27)，在t=tn和t=tn+1處取平均，且利用泰勒展開，可得到：

(28)

(29)

(30)

得到以下四階擬緊差分格式：

(31)

(32)

為了逼近Riesz空間分布階的分?jǐn)?shù)階擴散方程，將方程(31)改成：

(33)

線性系統(tǒng)(33)的系數(shù)矩陣A=(as,t)寫成下式：

(34)

則有以下的定理1。

定理1 差分格式(31)～(32)是唯一可解的。

3 穩(wěn)定性和收斂性的證明

為了證明差分格式的穩(wěn)定性和收斂性，有以下引理。

引理1A是自相似的，即對任意的u,v∈Vh,成立(Au,v)=(u,Av)。

證明

(35)

即證得。

證明

(36)

且

(37)

可得

(38)

即證得。

引理3 對任意的v∈Vh, 成立(δxv,v)≤0。

證明

(39)

使用引理3，容易得到以下定理2，其證明類似文獻(xiàn)[18]的定理2。

(40)

(41)

則有

(42)

由定理2可以直接推出如下定理3。

定理3 差分格式(31)～(32)按初值和右端項f是無條件穩(wěn)定的。

現(xiàn)在考慮差分格式 (31)～(32)的收斂性。

定理4 差分格式(31)～(32)是收斂的，且收斂階為O(τ2+σ4+h4)。

(43)

將式(28)和式(29)分別減去式(31)和式(32)，可得到誤差方程

(44)

(45)

根據(jù)定理2得到：

(46)

證畢。

4 數(shù)值例子

例1 考慮以下Riesz空間分布階擴散方程：

(47)

邊界條件：

u(0,t)=0,u(1,t)=0,t∈(0,1],

(48)

初始條件：

u(x,0)=x4(1-x)4,x∈[0,1]。

(49)

這里，

P(α)=-2Γ(9-α)cos(πα/2)，

(50)

且

(51)

方程(47)～(49)的精確解為u(x,t)=e-tx4(1-x)4。

圖1顯示了T=1.0時刻數(shù)值解與精確解的比較,取σ=1/1 000、τ=h2=1/400二者非常吻合。表 1 顯示σ=1/1 000、τ=h2、T=1.0時數(shù)值解與精確解的最大誤差及收斂。誤差率接近:Rε=lg(ε(h1)/ε(h2))/lg 2≈4，這與定理4中差分格式的收斂階是O(τ2+σ4+h4)是一致的。從圖1和表1可以看出，數(shù)值結(jié)果與理論分析結(jié)果是一致的。

表1 T=1.0時刻數(shù)值方法的最大誤差和收斂階(σ=1/1 000，τ=h2)

5 結(jié)論

本文發(fā)展了一個Riesz空間分布階的分?jǐn)?shù)階擴散方程的數(shù)值方法。將分布階方程離散為含有多項分?jǐn)?shù)階導(dǎo)數(shù)的微分方程，利用有限差分法對得到的微分方程進(jìn)行數(shù)值求解。證明差分格式是穩(wěn)定的和無條件收斂的。此外，本文的數(shù)值方法的收斂階為O(τ2+σ4+h4)。最后，給出了數(shù)值例子來證明本文數(shù)值方法與理論結(jié)果是一致的。這種方法和分析技術(shù)可用于求解和分析其他類型的分?jǐn)?shù)階偏微分方程。