微分求積方法解時(shí)間分?jǐn)?shù)階擴(kuò)散方程

曹 煥, 張學(xué)瑩, 劉 薈

(河海大學(xué) 理學(xué)院， 江蘇 南京 211100)

微分求積方法解時(shí)間分?jǐn)?shù)階擴(kuò)散方程

曹 煥, 張學(xué)瑩*, 劉 薈

(河海大學(xué) 理學(xué)院， 江蘇 南京 211100)

微分求積(DQ)法是一種基于徑向基函數(shù)(RBFs)插值的無網(wǎng)格方法.本文選MultiQ-uadrics(MQ)函數(shù)作為徑向基函數(shù),并采用微分求積方法解決時(shí)間分?jǐn)?shù)階擴(kuò)散方程。在離散過程中,采用有限差分法離散時(shí)間項(xiàng),采用微分求積方法離散空間項(xiàng).最后，結(jié)合數(shù)值求解的結(jié)果做出相應(yīng)的誤差分析.

徑向基函數(shù)； MQ函數(shù)； 微分求積法； 時(shí)間分?jǐn)?shù)階擴(kuò)散方程

0 引言

近年來,反常擴(kuò)散現(xiàn)象在物理學(xué)和工程學(xué)等領(lǐng)域得到了廣泛的關(guān)注[1,2],比如：污染物的運(yùn)輸、滲流、損耗以及磁等離子體等問題.與正常的擴(kuò)散現(xiàn)象相比,反常擴(kuò)散現(xiàn)象具有較強(qiáng)的遠(yuǎn)程相互作用與歷史依賴性等特點(diǎn),而標(biāo)準(zhǔn)的整數(shù)階偏微分方程已經(jīng)無法準(zhǔn)確的描述這類反常擴(kuò)散行為[3].隨后,分

數(shù)階導(dǎo)數(shù)被陸續(xù)的證實(shí)是一種可以準(zhǔn)確的描述這類反常擴(kuò)散現(xiàn)象的數(shù)學(xué)模型[4].所以,為了解決這些問題,引入了分?jǐn)?shù)階擴(kuò)散方程.

由于分?jǐn)?shù)階偏微分方程中含有分?jǐn)?shù)階的導(dǎo)數(shù),不容易求解分?jǐn)?shù)階偏微分方程的解析解,因此很多人研究其數(shù)值解法.目前,有限差分方法[5]已被證實(shí)為對(duì)分?jǐn)?shù)階擴(kuò)散方程進(jìn)行時(shí)間和空間上離散的一種有效方法.對(duì)于分?jǐn)?shù)階擴(kuò)散方程的數(shù)值解法，也可以選擇其他合適的方法對(duì)其空間上的離散,例如:邊界元法[6],有限元法[7,8],傅里葉法[9],無網(wǎng)格徑向基函數(shù)法[10-13]等.與傳統(tǒng)的有限差分方法對(duì)空間的離散相比較，對(duì)于處理一些大規(guī)模問題，這些方法在一定程度上減少了計(jì)算量，同時(shí)又提高了計(jì)算效率.2003年，Shu C提出DQ法解整數(shù)階的偏微分方程,且取得較理想的結(jié)果.本文將這一方法運(yùn)用到解分?jǐn)?shù)階偏微分方程中，并給出誤差結(jié)果分析.

1 基本理論

1.1 時(shí)間分?jǐn)?shù)階擴(kuò)散方程

(1)

0lt;αlt;1,X∈Ω,t∈(0,T)

邊界條件:

u(X,t)=g(X,t)，X∈?Ω,t∈(0,T)

(2)

初始條件:

u(X,0)=u0(X),X∈Ω

(3)

式(1)～(3)中：Q(X,t),g(X,t),u0(x)是給定的函數(shù),T是總時(shí)間,?α/?tα是關(guān)于時(shí)間t的分?jǐn)?shù)階導(dǎo).

本文采用Caputo時(shí)間分?jǐn)?shù)階導(dǎo)數(shù)的定義,定義如下:

(4)

這里0lt;αlt;1.

1.2 FDM方法對(duì)時(shí)間的離散

首先把區(qū)間[0,T]分割成K個(gè)小區(qū)間,定義tk=kτ,k=0,1,…,K,這里τ=T/K為時(shí)間步長(zhǎng).用u(X,tk)表示函數(shù)u(X,t)在tk處的精確解.本文利用有限差分方法在t=tk+1處對(duì)時(shí)間分?jǐn)?shù)階導(dǎo)數(shù)進(jìn)行離散,有

(5)

式(5)中：

根據(jù)式(5),對(duì)方程式(1)進(jìn)行時(shí)間離散,得到

Δuk+1+Qk+1=

(6)

1.3 微分求積方法對(duì)空間的離散

設(shè)f(x)是光滑函數(shù),其中x=(x1,x2,…,xn)T.已知它在節(jié)點(diǎn)xi及其支撐域內(nèi)各支撐點(diǎn)xj,j=1,2,…,n上的函數(shù)值,則f(x)關(guān)于xk的m階導(dǎo)數(shù)可以表示為

(7)

(8)

可將其寫成矩陣形式是

(9)

由文獻(xiàn)[14]可知式(9)中的矩陣是條件正定的,故系數(shù)矩陣是可逆的,從而可求得系數(shù).再將其代入式(7),就可以求出函數(shù)f(x)關(guān)于xk的m階導(dǎo)數(shù)的近似值.

由上述過程,本文采用DQ方法對(duì)方程式(6)進(jìn)行空間離散,可知

(10)

再將式(10)代入式(6)得到對(duì)空間的離散,并對(duì)其整理可得

(11)

2 數(shù)值算例

本文考慮在規(guī)則區(qū)域Ω=[0,0.5]×[0,0.5]內(nèi),檢驗(yàn)時(shí)間分?jǐn)?shù)階擴(kuò)散方程的有效性和準(zhǔn)確性.數(shù)值的準(zhǔn)確性由均方根誤差(RMSE),最大絕對(duì)誤差(MAE)和最大相對(duì)誤差(RAE)計(jì)算.其誤差定義如下：

算例

0lt;αlt;1,X∈Ω,t∈T

邊界條件:

u(X,t)=t2ex+y,

X∈?Ω,t∈(0,T)

初始條件:

u(X,0)=0,X∈Ω.

α=0.9,Ω=[0,0.5]2，

Q(X,t)=(2t2-α/Γ(3-α)-2t2)ex+y.

該方程的精確解為

u(X,t)=t2ex+y,x∈Ω,t∈(0,T).

圖1為選取網(wǎng)格點(diǎn)N=961,分?jǐn)?shù)階導(dǎo)數(shù)α=0.9,時(shí)間步長(zhǎng)dt=0.1,參數(shù)c=0.05時(shí)的絕對(duì)誤差和相對(duì)誤差圖.

(a) N=961,α=0.9,dt=0.1,c=0.5時(shí)絕對(duì)誤差圖

(b) N=961,α=0.9,dt=0.1,c=0.5時(shí)相對(duì)誤差分布圖

圖2為選取分?jǐn)?shù)階導(dǎo)數(shù)α=0.9,參數(shù)c=0.02時(shí),隨著插值點(diǎn)個(gè)數(shù)增加,不同時(shí)間步長(zhǎng)下的均方差對(duì)比圖;圖3為選取分?jǐn)?shù)階導(dǎo)數(shù)α=0.9,時(shí)間步長(zhǎng)t=0.1時(shí),隨著插值點(diǎn)個(gè)數(shù)增加,不同參數(shù)下的均方差對(duì)比圖.

圖2 不同時(shí)間步長(zhǎng)下的均方根誤差比較

圖3 不同參數(shù)c下的均方根誤差比較

從圖2和圖3可以看出,參數(shù)c和時(shí)間步長(zhǎng)的選取對(duì)數(shù)值結(jié)果有較大的影響,當(dāng)選取較大的時(shí)間步長(zhǎng)時(shí),精度有明顯的提高,當(dāng)插值點(diǎn)個(gè)數(shù)為256時(shí),精度達(dá)到最高,當(dāng)插值點(diǎn)個(gè)數(shù)的增加,精度變化幅度逐漸減小,最終趨于穩(wěn)定,但如果點(diǎn)再繼續(xù)增加,就可能會(huì)出現(xiàn)病態(tài);當(dāng)選取較小的參數(shù)c時(shí),隨著插值點(diǎn)個(gè)數(shù)的增加,精度開始有明顯的提高,但當(dāng)插值點(diǎn)個(gè)數(shù)再增加時(shí),精度變化幅度逐漸減小,最終趨于穩(wěn)定.因此,選取適當(dāng)?shù)臅r(shí)間步長(zhǎng)和參數(shù)c是提高數(shù)值結(jié)果精度的有效方法.

表1為選取時(shí)間步長(zhǎng)t=0.1,參數(shù)c=0.02時(shí),不同分?jǐn)?shù)階導(dǎo)數(shù)下的絕對(duì)誤差和相對(duì)誤差的結(jié)果比較.

表1 t=0.1，c=0.02時(shí)不同α下的絕對(duì)誤差和相對(duì)誤差比較

從表1可以看出,當(dāng)取時(shí)間步長(zhǎng)t=0.1,參數(shù)c=0.02,分?jǐn)?shù)階的導(dǎo)數(shù)越接近1時(shí),隨著插值點(diǎn)個(gè)數(shù)的增加,誤差越來越精確.

3 結(jié)論

本文在結(jié)合有限差分方法的基礎(chǔ)上，應(yīng)用微分求積方法對(duì)規(guī)則區(qū)域上的時(shí)間分?jǐn)?shù)階擴(kuò)散方程進(jìn)行數(shù)值求解.根據(jù)Caputo分?jǐn)?shù)階導(dǎo)數(shù)的積分定義,利用有限差分方法對(duì)時(shí)間導(dǎo)數(shù)進(jìn)行離散,利用微分求積方法對(duì)空間導(dǎo)數(shù)進(jìn)行離散,并通過數(shù)值算例得出數(shù)值結(jié)果.

數(shù)值結(jié)果表明,微分求積方法是解決時(shí)間分?jǐn)?shù)階擴(kuò)散方程比較可行的方法.在解決問題的過程中，選取了MQ函數(shù)作為徑向基函數(shù),從圖3可以看出參數(shù)c的選擇對(duì)數(shù)值結(jié)果具有很大的影響,當(dāng)選取適當(dāng)參數(shù)c時(shí),數(shù)值結(jié)果的精度更高.同時(shí)，還可得出時(shí)間步長(zhǎng)和插值點(diǎn)個(gè)數(shù)也對(duì)其具有一定影響.

后續(xù)，還可進(jìn)一步討論在不規(guī)則區(qū)域上此方法的可行性以及怎樣選取參數(shù)c可使數(shù)值結(jié)果達(dá)到更好的精度.

[1] R Gorenflo,F Mainardi,D Noretti,et al.Discrete random walk models for space time fractional diffusion[J].Chenmical Physics,2007,284(1-2):521-541.

[2] L Baowen,W Jiao.Anomalous heat conduction and anomalous diffusion in onedimensional systems[J].Physical Review Letters,2003,91(4):318-322.

[3] Zhang X Y,Zhu H Q,Kuo L H.A comparison study of the LMAPS method and the LDQ method for time-dependent problems[J].Engineering Analysis with Boundary Elements,2013,37:1 408-1 415.

[4] F Liu,V Anh,I Turner,et al.Time fractional advection dispersion equation[J].Journal of Applied Mathematics and Computing,2003,13:233-245.

[5] H Wang,K Wang,T Sircar.A direct O (Nlog2N) finite difference method for fractional diffusion equations[J].Journal of Computational Physics,2010,229(299):8 095-8 104.

[6] Katsikadelis J T.The BEM for numerical solution of partial fractional differential equations[J].Computers amp; Mathematics with Applications,2011,62(3):891-901.

[7] Li C,Zhao Z, Chen Y. Numerical approximation of nonlinear fractional differential equations with subdiffusion and superdiffusion[J].Computers amp; Mathematics with Applications,2011,62(3):855-875.

[8] Jiang Y,Ma J.High-order finite element methods for time-fractional partial differential equations[J].Journal of Computational amp; Applied Mathematics,2011,235(11):3 285-3 290.

[9] Chen C M,Liu F,Turner I,et al.A Fourier method for the fractional diffusion equation describing sub-diffusion[J].Journal of Computational Physics,2007,227(2):886-897.

[10] Chen W,Ye L,Sun H.Fractional diffusion equations by the Kansa method[J].Computers amp; Mathematics with Applications,2010,59(5):1 614-1 620.

[11] Liu Q,Gu Y,Zhang P,et al.An implicit RBF meshless approach for time fractional diffusion equations[J].Computational Mechanics,2011,48(1):1-12.

[12] Fu Z J,Wen Chen,Leevan Ling,el al.Method of approximate particular solutions for constantand variable-order fractional diffusion modes[J].Engineering Analysis with Boundary Elements,2015,57:37-46.

[13] Chen W,Fu Z J.Recent advances on radial basis function collocation methods[M].Berlin:Springer Berlin Heidelberg,2014.

[14] Micchelli C A.Interolation of scattered deta:Distance matrices and conditionally positive definite functions[J].Constr Approx,1986(2):11-22.

【責(zé)任編輯：蔣亞儒】

Thedifferentialquadraturemethodsolutionsforthetime-fractionaldiffusionequation

CAO Huan, ZHANG Xue-ying*, LIU Hui

(School of Science, Hohai University, Nanjing 211100, China)

The differential quadrature method is an alternative radial basis functions (RBFs) meshless method.This article selects MQ function as the radial basis function,and applys the differential quadrature method to solve the time fractional diffusion equation.In the discretization formulation,a finite difference scheme and the DQ are used repectively to discretize time fractional derivative and spatial derivative terms.Finally we make the error analysis with the results of the numerical investigation example.

RBFs; MQ function; DQ method; time fractional diffusion function

2017-08-27

教育部留學(xué)回國(guó)人員科研啟動(dòng)基金項(xiàng)目(20145003412)； 江蘇省自然科學(xué)基金項(xiàng)目(BK20160853)

曹 煥(1992-)，女，山東菏澤人，在讀碩士研究生，研究方向：計(jì)算數(shù)學(xué)

張學(xué)瑩(1973-)，男，山東郯城人，副教授，碩士生導(dǎo)師，研究方向：計(jì)算數(shù)學(xué)， zhangxy@hhu.edu.cn

2096-398X(2017)06-0179-04

O241.82

A