吳 迪， 李小林

（重慶師范大學(xué) 數(shù)學(xué)科學(xué)學(xué)院，重慶 401331）

引 言

時(shí)間分?jǐn)?shù)階擴(kuò)散波方程可以用來(lái)模擬一些重要的物理現(xiàn)象[1-2]，比如黏彈性介質(zhì)中機(jī)械波的傳播、具有記憶的非Markov 擴(kuò)散過(guò)程、以及非晶態(tài)半導(dǎo)體中的電荷運(yùn)輸.有限差分法[1-5]和有限元法[5-8]等方法在近年來(lái)被用于數(shù)值求解此方程.

無(wú)網(wǎng)格法基于點(diǎn)的近似，擺脫了有限差分和有限元等方法中網(wǎng)格單元的限制，具有實(shí)施便利和精度高等優(yōu)勢(shì)[9].Yang 等[10]和Kumar 等[11]發(fā)展了數(shù)值求解時(shí)間分?jǐn)?shù)階擴(kuò)散波方程的無(wú)網(wǎng)格配點(diǎn)法.無(wú)單元Galerkin法[8,12-13]是一種基于移動(dòng)最小二乘近似的全局弱式無(wú)網(wǎng)格法，比無(wú)網(wǎng)格配點(diǎn)法具有更好的穩(wěn)定性.Dehghan 和Abbaszadeh 通過(guò)使用L1 逼近公式離散時(shí)間變量，在文獻(xiàn)[14]中建立了數(shù)值求解時(shí)間分?jǐn)?shù)階擴(kuò)散波方程的無(wú)單元Galerkin 法，但是只處理了Neumann 邊界條件，沒有包含無(wú)單元Galerkin 法不易處理的Dirichlet 邊界條件.他們?cè)谖墨I(xiàn)[8]和[15]中用插值型無(wú)單元Galerkin 法，數(shù)值求解了含有Dirichlet 邊界條件的時(shí)間分?jǐn)?shù)階擴(kuò)散波方程，但是在構(gòu)造形函數(shù)時(shí)需要使用特殊的奇異權(quán)函數(shù).

本文研究了Caputo 意義下的時(shí)間分?jǐn)?shù)階擴(kuò)散波方程的無(wú)單元Galerkin 法.首先，類似于文獻(xiàn)[8,14-15]，基于文獻(xiàn)[1]中的L1 逼近公式，離散時(shí)間變量.其次，為了有效處理Dirichlet 邊界條件并避免使用奇異權(quán)函數(shù)，本文通過(guò)罰函數(shù)方法處理Dirichlet 邊界條件，建立了該方程的無(wú)單元Galerkin 法離散線性代數(shù)方程組.然后，借鑒文獻(xiàn)[16-17]的技巧理論分析了時(shí)間分?jǐn)?shù)階擴(kuò)散波方程的無(wú)單元Galerkin 法的誤差.最后，利用數(shù)值算例，進(jìn)一步論證了該方法的有效性和精度.

1 問題描述

考慮如下時(shí)間分?jǐn)?shù)階擴(kuò)散波方程[1]：

初始條件為

邊界條件為

其中c為擴(kuò)散波動(dòng)系數(shù)， Δu=?2u/?x2+?2u/?y2， Ω 是問題區(qū)域， Γ1是 Dirichlet 邊界， Γ2是 Neumann 邊界，n是邊界?Ω=Γ1∪Γ2上 的外法向單位向量，f(x,y)， ψ (x,y) ，φ (x,y)， ξ (x,y)， γ (x,y)是已知函數(shù)，且

是Caputo 意義下的分?jǐn)?shù)階導(dǎo)數(shù)，Γ (·)為Gamma 函數(shù).

2 無(wú)單元Galerkin 法數(shù)值離散

2.1 時(shí)間離散

令τ 表示時(shí)間步長(zhǎng)，un=u(x,y,nτ)，則Caputo 分?jǐn)?shù)階導(dǎo)數(shù)可利用L1 逼近公式[1, 3]近似為

將?αun/?tα和?αun?1/?tα的表達(dá)式代入式(7)，并利用式(6)可得

2.2 Galerkin 變分方程

橢圓邊值問題(10)可用無(wú)單元Galerkin 法迭代求解.由于移動(dòng)最小二乘近似生成的形函數(shù)缺乏插值性質(zhì)，用無(wú)單元Galerkin 法求解問題(10)時(shí)，Dirichlet 邊界條件的施加較困難[18]，可以采用罰函數(shù)法來(lái)處理.引入罰因子β，則橢圓邊值問題(10)可近似為

2.3 無(wú)網(wǎng)格離散

這里

其中 M是移動(dòng)最小二乘近似中的逼近算子.為了便于數(shù)值計(jì)算，由Gauss 公式和式(11)中的邊界條件，進(jìn)一步可得

將式(18)代入式(16)，可得

由v的任意性，在式(19)中讓v依 次取 φ1， φ2， ···， φN，則橢圓邊值問題(10)的無(wú)網(wǎng)格離散為

其中

同理可得當(dāng)n=1時(shí)，橢圓邊值問題(10)可離散為如下線性代數(shù)方程組：

其中

由方程組(21)可以依次迭代求解得到Un.最終，由式(17)可以求出時(shí)間分?jǐn)?shù)階擴(kuò)散波方程的無(wú)單元Galerkin 法近似解.

3 誤差分析

定理1設(shè)un∈Hq+1(Ω) 是 問題(10)的解析解，為近似變分問題(16)的數(shù)值解，則時(shí)間分?jǐn)?shù)階擴(kuò)散波方程的無(wú)單元Galerkin 法有如下誤差：

所以

由文獻(xiàn)[16]可得

4 數(shù)值算例

4.1 算例1

考慮如下時(shí)間分?jǐn)?shù)階擴(kuò)散波方程[8]：

初始條件和邊界條件分別為

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

u(x,y,t)=t2sin(x+y), (x,y)∈?Ω,t∈[0,1].

該問題的解析解為u(x,y,t)=t2sin(x+y).

圖1 給出了 α =1．65，h=1/50和 τ =1/40 時(shí)，無(wú)單元Galerkin 法獲得的近似解Uh和 誤差 |Uh?u|，數(shù)據(jù)顯示誤差小于 1 ．5×10?4，說(shuō)明該方法的計(jì)算精度較高.為了研究罰因子 β 對(duì)誤差的影響，圖2 給出了 α=1．65，τ=1/200 時(shí) ，常數(shù)罰因子 β 和變化罰因子 β =Cβh?2對(duì)誤差和收斂性的影響.罰因子β 取為常數(shù)時(shí)，太大或太小都會(huì)讓誤差增加；罰因子 β =Cβh?2隨著空間步長(zhǎng)h變 化時(shí)，誤差始終隨著空間步長(zhǎng)h的減小而減小，在后面的計(jì)算中選取罰因子β =102h?2.

圖1 算例1 在α =1．65，h =1/50 和τ =1/40時(shí)的近似解和誤差：(a) 近似解；(b) 誤差Fig.1 Graphs of approximate solutions and resulting errors with α = 1．65, h =1/50 and τ =1/40 in example 1: (a) the approximate solutions;(b) the resulting errors

圖2 常數(shù)罰因子β 和變化罰因子β =Cβh?2 對(duì)誤差和收斂性的影響：(a) 常數(shù)罰因子β ；(b) 變化罰因子β=Cβh?2Fig.2 The error and convergence for fixed penalty factor β and variable penalty factor β =Cβh?2 : (a) for fixed penalty factor β ;(b) for variable penalty factor β=Cβh?2

當(dāng) α =1．65， 1.75，1.85，圖3(a)給出了h=1/20 時(shí) ，無(wú)單元Galerkin 法關(guān)于時(shí)間步長(zhǎng)τ 和 空間步長(zhǎng)h的收斂性；圖3(b)給出了 τ=1/1 000時(shí) ，無(wú)單元Galerkin 法關(guān)于空間步長(zhǎng)h的 收斂性.圖3(a)使用的時(shí)間步長(zhǎng) τ為1/8，1/16，1/32，1/64；圖3(b)使用的空間步長(zhǎng)h為1/2，1/4，1/8，1/16，對(duì)應(yīng)的節(jié)點(diǎn)個(gè)數(shù)為9，25，81，289.圖3(a)和(b)中的數(shù)值結(jié)果分別表明，無(wú)單元Galerkin 法關(guān)于時(shí)間步長(zhǎng)和空間步長(zhǎng)都得到了很好的收斂結(jié)果.

圖3 算例1 關(guān)于時(shí)間步長(zhǎng)τ 和空間步長(zhǎng)h 的收斂性：(a) 時(shí)間步長(zhǎng)τ；(b) 空間步長(zhǎng)hFig.3 The convergence with respect to time step τ and spatial step h in example 1: (a) for time step τ; (b) for spatial step h

影響域半徑可設(shè)置為Sscale·h[9,17]，其中Sscale是控制影響域大小的參數(shù).當(dāng) τ =1/10 和h=1/20時(shí)，圖4 展示了Sscale對(duì)誤差的影響，結(jié)果表明，影響域大小對(duì)誤差的影響較小，在后面的計(jì)算中選取Sscale=2.5.表1 比較了無(wú)單元Galerkin 法(EFG)和有限元法(FEM)[8]在不同分?jǐn)?shù)階情況下的L∞誤差，結(jié)果表明無(wú)單元Galerkin 法具有更高的計(jì)算精度.

圖4 影響域參數(shù)Sscale 對(duì)誤差的影響Fig.4 Effects of influence domain parameter Sscale on errors

表1 比較有限元法和無(wú)單元Galerkin 法的 L∞誤差Table 1 Comparison of L∞ errors between the finite element and the element-free Galerkin methods

4.2 算例2

考慮如下時(shí)間分?jǐn)?shù)階擴(kuò)散波方程[8]：

初始條件和邊界條件分別為

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

u(x,y,t)=0, (x,y)∈?Ω,t∈[0,1].

該問題的解析解為u(x,y,t)=t2+αsin(x)sin(y).

當(dāng) α =1．5時(shí) ，圖5(a)和(b)分別給出了h=π/40和 τ =1/256， 無(wú)單元Galerkin 法(EFG)關(guān)于時(shí)間步長(zhǎng)τ 和空間步長(zhǎng)h的收斂性.為了比較，圖5 同時(shí)給出了Galerkin 有限元法(Galerkin FEM)[8]和交替方向隱式有限元法(ADI FEM)[6]的計(jì)算結(jié)果.數(shù)值結(jié)果表明無(wú)單元Galerkin 法得到了很好的收斂結(jié)果，并且比兩種有限元法的計(jì)算精度更高.表2 中比較了當(dāng) α=1．25，h=π/40時(shí)，Galerkin 有限元法[8]、交替方向隱式有限元法[6]和無(wú)單元Galerkin 法的L2誤差，顯然無(wú)單元Galerkin 法計(jì)算精度更高.

圖5 算例2 關(guān)于時(shí)間步長(zhǎng)τ 和空間步長(zhǎng)h 的收斂性：(a) 時(shí)間步長(zhǎng)τ；(b) 空間步長(zhǎng)hFig.5 The convergence with respect to time step τ and spatial step h in example 2: (a) for time step τ; (b) for spatial step h

表2 比較h =π/40時(shí) ，三種方法的 L2誤差Table 2 Comparison of L2 errors obtained from 3 numerical methods with h=π/40

5 結(jié) 論

本文構(gòu)造了數(shù)值求解時(shí)間分?jǐn)?shù)階擴(kuò)散波動(dòng)方程的無(wú)單元Galerkin 法，詳細(xì)推導(dǎo)了該方法的計(jì)算公式和理論誤差.誤差理論表明，數(shù)值解的誤差與時(shí)間步長(zhǎng)τ、空間步長(zhǎng)h和 罰因子β 有關(guān)，并且在時(shí)間方向和空間方向的收斂速率分別約為O(τ3?α)和O(h2).數(shù)值算例驗(yàn)證了該方法的有效性和理論分析結(jié)果，并且數(shù)值結(jié)果表明本文方法比有限元等方法擁有更高的計(jì)算精度.本文采用了L1 公式對(duì)時(shí)間變量進(jìn)行離散，未來(lái)還可以研究利用L(1-2)公式、加權(quán)位移的G-L 公式等高階數(shù)值解法去離散時(shí)間變量，同時(shí)還可以嘗試結(jié)合更多的無(wú)網(wǎng)格方法去構(gòu)建求解時(shí)間分?jǐn)?shù)階擴(kuò)散波方程的時(shí)空高階格式.另一方面，本文所研究的時(shí)間分?jǐn)?shù)階僅為常數(shù)，未來(lái)還可以對(duì)變時(shí)間分?jǐn)?shù)階微分方程的數(shù)值求解進(jìn)行研究.