朱瑞,張根根,肖飛雁,蘭海峰

(廣西師范大學(xué)數(shù)學(xué)與統(tǒng)計(jì)學(xué)院,廣西 桂林541004)

1.引言

分?jǐn)?shù)階延遲微分方程(FDDEs)在計(jì)算機(jī)神經(jīng)網(wǎng)絡(luò)[19],流體力學(xué)[18],電磁學(xué)[11],工程學(xué)[20?21]等應(yīng)用科學(xué)領(lǐng)域得到廣泛應(yīng)用.關(guān)于分?jǐn)?shù)階延遲微分方程定性理論,Benchohra等[1]首先討論了一類帶無窮時(shí)滯的Riemann-Liouville分?jǐn)?shù)階微分方程初值問題解的存在性.DENG[2]研究了一類分?jǐn)?shù)階延遲微分方程解的存在唯一性并得到兩個(gè)新的存在性和唯一性結(jié)果.在分?jǐn)?shù)階延遲微分方程數(shù)值求解方面,Moghaddam[3]提出了利用分?jǐn)?shù)階有限差分方法求解分?jǐn)?shù)階延遲微分方程.Mogrado等[4]利用自適應(yīng)的分?jǐn)?shù)階向后差分方法求解了一類常延遲線性分?jǐn)?shù)階微分方程初值問題.Bhalekar等[5]將Adams-Bashforth-Moulton算法擴(kuò)展到求解分?jǐn)?shù)階延遲微分方程.WANG[6]又將Adams-Bashforth-Moulton方法與線性插值方法結(jié)合來近似分?jǐn)?shù)階延遲微分方程.Daftardargejji[7]提出了一種新的用于分?jǐn)?shù)階延遲微分方程的預(yù)估校正方法.

另一方面,關(guān)于卷積Runge-Kutta方法,Lubich[8?9,13?14]將卷積核進(jìn)行Lapace變換,利用Runge-Kutta 方法構(gòu)造了求解分?jǐn)?shù)階微分方程的數(shù)值算法,并給出相應(yīng)的理論分析結(jié)果.接著,在文[10]中給出了基于A-穩(wěn)定Runge-Kutta方法的卷積積分的誤差分析,并證明該方法近似的階數(shù)取決于Runge-Kutta方法的經(jīng)典階數(shù)和Lapace變換的增長(zhǎng)指數(shù).曹學(xué)年等[12]構(gòu)造了求解非線性分?jǐn)?shù)階微分方程的Radau IIA方法,并證明該方法的相容性,收斂性和穩(wěn)定性.徐大[22]利用Runge-Kutta方法研究Volterra型積分方程初邊值問題,并證明了數(shù)值方法的收斂性.

本文基于強(qiáng)A-穩(wěn)定Runge-Kutta方法,構(gòu)造了求解非線性分?jǐn)?shù)階延遲微分方程的離散格式.第三節(jié)和第四節(jié)分別對(duì)數(shù)值方法進(jìn)行了誤差分析和穩(wěn)定性分析.最后給出兩個(gè)數(shù)值算例,以驗(yàn)證提出方法的有效性.

2.數(shù)值格式

考慮如下分?jǐn)?shù)階延遲微分方程初值問題

其中常延遲τ ≥0,φ(k)∈C[?τ,0]為初值函數(shù),f:[0,T]×Rd×Rd→Rd是一個(gè)連續(xù)映射且滿足經(jīng)典Lipschiz條件

這里L(fēng)1,L2為常數(shù),符號(hào)∥·∥表示空間Rd中的范數(shù).

定義2.1[15]表示Caputo分?jǐn)?shù)階導(dǎo)數(shù),其定義形式為:

其中l(wèi) ?1<α ≤l,l ∈N,t>0.

定義2.2[15]表示Riemann-Liouville分?jǐn)?shù)階導(dǎo)數(shù),其定義形式為:

根據(jù)Riemann-Liouville分?jǐn)?shù)階導(dǎo)數(shù)與Caputo分?jǐn)?shù)階導(dǎo)數(shù)的關(guān)系,微分方程(2.1)等價(jià)于以下方程

并且記問題(2.3)的真解為y(t).

應(yīng)用基于Riemann-Liouville分?jǐn)?shù)階導(dǎo)數(shù)的強(qiáng)A-穩(wěn)定Runge-Kutta方法[9]

求解(2.3)式,得到在點(diǎn)tn+1=(n+1)h處的離散格式

即

其中

3.誤差分析

在給出主要結(jié)果之前,我們首先介紹幾個(gè)有用的引理.

引理3.1[9]假設(shè)K(s)在扇形區(qū)域|arg(s ?c)| <π??,其中? <π2,c ∈R,存在實(shí)數(shù)μ,M,使得|K(s)|≤M ·|s|?μ,考慮一個(gè)強(qiáng)A穩(wěn)定的m級(jí)p階Runge-Kutta方法,此時(shí)存在一個(gè)h0>0,當(dāng)0

其中,q為級(jí)階,特別地,對(duì)于Radau IIA方法有p=2m ?1,q=m.

引理3.2[9]在引理3.1的假設(shè)條件下,存在h0>0,使得當(dāng)h

其中n=0與n=1有相同的界限,c,γ,h0僅依賴于方法和引理(3.1)中假設(shè)條件中的常數(shù),可以獨(dú)立于μ.

引理3.3(Gronwall不等式) 令a,b ≥0,h>0,ηi滿足

則必有

其中M0=max(|η0|,|η1|,··· ,|ηk?1|).

接下來,我們將給出Runge-Kutta方法的局部截?cái)嗾`差和收斂性分析.記Un為mN維向量,U(1)n=O(hm+1)為Un前(m ?1)N個(gè)分量,U(2)n=O(hp)+O(hm+1+μ|logh|)為后N個(gè)分量.

定理3.1對(duì)于分?jǐn)?shù)階延遲微分方程(2.3),基于m級(jí)Runge-Kutta方法的局部截?cái)嗾`差為O(hp)+O(hm+1+μ|logh|)).

證假設(shè)Yvj=y(tv+cjh)(v ≤n).由(2.3),(2.6) 及引理3.1 可得

這里ω0,m是矩陣W0=h?α?1(A?1)α的第m行的第m個(gè)元素,令hω0,m=h?α,其中為依賴于方法的非零常數(shù).將hα/乘以上式,兩邊取范數(shù)可推出

存在最大步長(zhǎng)h0>0,使得則當(dāng)0

證畢.

定理3.2對(duì)于分?jǐn)?shù)階延遲微分方程(2.3),基于強(qiáng)A-穩(wěn)定的m級(jí)Runge-Kutta方法滿足

證記evj=y(tv+cjh)?Yvj,j=1,··· ,m,Ev=(ev1,ev2,··· ,evm)T,由(2.3),(2.6)及引理3.1可得

即

其中,W0=h?α?1(A?1)α為非奇異矩陣.將上式左右兩端同時(shí)左乘hαAα有

注意到

對(duì)于∥En?u∥進(jìn)行如下討論.

情況I當(dāng)n?u ≤0時(shí)，y(t?τ)=φ(t?τ)根據(jù)(2.4)式可知0.對(duì)(3.2)式兩邊取范數(shù),有

存在h1>0使得(hα1∥Aα∥)?1>L1,則當(dāng)0

這里C1=max(∥Wv∥,v=0,··· ,n ?1),

由引理3.3有

當(dāng)m ≥2時(shí),∥Un∥=O(hm+1)；當(dāng)m=1時(shí),∥Un∥=O(hp),所以

情況II當(dāng)n ?u ≥0時(shí),

證畢.

4.穩(wěn)定性分析

定理4.1卷積Runge-Kutta方法求解分?jǐn)?shù)階延遲微分方程(2.3)是穩(wěn)定的.

證設(shè){yn+1},{xn+1}分別是由分?jǐn)?shù)階延遲微分方程(2.3)及其擾動(dòng)問題

的解,令

其中xn+1=Xnm,yn+1=Ynm.

記σij=Xij?Yij,j=1,··· ,m,?v=(σv1,···,σvm)T,v=?u,?u+1,··· ,n,?k(0)=xk(0)?yk(0),k=0,1,··· ,l ?1,則

即

注意到

且

對(duì)(4.1)式兩邊同乘hαAα,并兩邊取范數(shù)得

即

取h2為最大步長(zhǎng),且使(hα2∥Aα∥)?1>L1,則當(dāng)0

即

這里C3,C4為常數(shù),根據(jù)引理3.3有

從而

所以數(shù)值方法穩(wěn)定.證畢.

5.數(shù)值試驗(yàn)

例5.1考慮如下分?jǐn)?shù)階延遲微分方程問題

其精確解為y=t3.

利用三階Radau IIA方法,分別取不同步時(shí)刻t和不同α在t ∈[0,2]進(jìn)行數(shù)值試驗(yàn),其絕對(duì)誤差e如表5.1所示.取h=1/20,α=0.1,0.3,0.5,0.7,0.9時(shí)的誤差分布圖見圖5.1.

表5.1 在t=2.0時(shí)刻,不同步長(zhǎng)h和α下數(shù)值解的絕對(duì)誤差

圖5.1 h=1/20時(shí)取不同α值的誤差分布

圖5.2 h=1/40時(shí)取不同α值的誤差分布

例5.2考慮如下分?jǐn)?shù)階延遲微分方程問題

其真解為y=t2.

利用三階Radau IIA方法,分別取不同步長(zhǎng)h和不同α在t ∈[0,2]進(jìn)行數(shù)值試驗(yàn),其絕對(duì)誤差如表5.2所示.取h=1/40,α=0.1,0.3,0.5,0.7,0.9時(shí)的誤差分布圖見圖5.2.

表5.2 在t=2.0時(shí)刻,不同步長(zhǎng)h和α下數(shù)值解的絕對(duì)誤差