求解分?jǐn)?shù)階常微分方程的一個(gè)高階數(shù)值逼近格式*

曹文平，肖承家，王自強(qiáng)

(貴州民族大學(xué) 數(shù)據(jù)科學(xué)與信息工程學(xué)院，貴州 貴陽(yáng) 550025)

0 引言

20世紀(jì)70年代末， Mandelbrot[1]提出了分形學(xué)說(shuō)，并且將Riemann-Liouville分?jǐn)?shù)階微積分用以分析和研究分形媒介中的布朗運(yùn)動(dòng)以后，分?jǐn)?shù)階算子尤其是分?jǐn)?shù)階微積分和分?jǐn)?shù)階微分方程理論及其應(yīng)用研究在國(guó)際上才得到廣泛關(guān)注和迅速發(fā)展。

與整數(shù)階微分方程不同，分?jǐn)?shù)階微分方程的理論研究在文獻(xiàn)中很少看到。Diethelm等人考慮了分?jǐn)?shù)階常微分方程(FODEs)初值問(wèn)題解的適定性[2]， Kilbas等人利用廣義的Mittag-Leffler函數(shù)研究了Volterra integro-differential方程的解的表達(dá)式[3]，Diethelm[4]給出了FODEs理論方面的最新發(fā)展。

給定一個(gè)一般的右端項(xiàng)f，確定分?jǐn)?shù)階微分方程的解析解是相當(dāng)困難的。因此，我們必須尋找有效求解FODEs的數(shù)值方法。

我們考慮分?jǐn)?shù)階常微分方程如下：

(1)

且滿足如下的初值條件：

y(0)=y0

(2)

(3)

其中表達(dá)式Γ(·)表示Euler Gamma函數(shù)。

對(duì)于方程(1)的數(shù)值方法，許多研究者們已經(jīng)做了大量的工作[6]。文獻(xiàn)[7，8]將(1)式轉(zhuǎn)化為Volterra型積分方程，則可利用積分方程的一些技巧建立了高階數(shù)值格式。本文中，我們利用二次拉格朗日插值來(lái)逼近分?jǐn)?shù)階導(dǎo)數(shù)，得到了一個(gè)高階數(shù)值逼近格式。

本文第一部分， 詳細(xì)地構(gòu)造了一個(gè)高階數(shù)值逼近格式; 第二部分，列出了所構(gòu)造數(shù)值格式的局部截?cái)嗾`差估計(jì); 第三部分，用一系列算例來(lái)驗(yàn)證理論預(yù)測(cè)的正確性。

1 高階數(shù)值格式

為了構(gòu)造高階數(shù)值逼近格式，我們將區(qū)間[0，T]分成2N個(gè)等分的子區(qū)間，設(shè)：

h=T/(2N)，記xj=jh，j=0，1，…，2N

在xj處的數(shù)值解記為yi，且

下面，我們開(kāi)始構(gòu)造高階數(shù)值逼近格式。

(4)

對(duì)于子區(qū)間[x2k，x2k+2]，k=0，1，…，m-1，y(x)的逼近形式如下：

(5)

其中φi，k(x)，i=0，1，2；k=0，1，…，m-1，是定義在點(diǎn)x2k，x2k+1，x2k+2上的二次拉格朗日基函數(shù)：

(6)

對(duì)于子區(qū)間[x2m，x2m+1]，y(x)的逼近形式為：

(7)

(8)

于是，

[I[x2m，x2m+1]y(s)]′ds

(9)

其中：

i=0，1，2；k=0，1，…m-1

(10)

i=0，1，2

可知，(10)式是可以精確計(jì)算出來(lái)的。y2m+1/2通過(guò)下面的插值來(lái)獲得：

(11)

因此，第2m+1步的數(shù)值格式如下：

(12)

(13)

對(duì)于[x2k，x2k+2]，k=0，1，…，m，y(x)的逼近式見(jiàn)(5)式，且φi，k(x)，i=0，1，2；k=0，1，…，m，見(jiàn)(6)式。于是，

[I[x2k，x2k+2]y(s)]′ds

其中：

因此我們得出2m+2步的數(shù)值格式:

(14)

綜上所述，結(jié)合(12)式和(14)式，我們得到數(shù)值格式如下：

(15)

2 局部截?cái)嗾`差估計(jì)

數(shù)值格式(15)的局部截?cái)嗾`差估計(jì)如下。

定理1.假定y(x)∈C4[0，T].設(shè)：

(16)

|rk(h)|≤Ch3-α，k=1，2，…，2N

(17)

其中C是一個(gè)依賴于y但與h無(wú)關(guān)的常數(shù)。

3 數(shù)值算例

我們要進(jìn)行一系列數(shù)值試驗(yàn)，來(lái)測(cè)試格式的有效性。準(zhǔn)確地說(shuō)，我們的主要目的是驗(yàn)證數(shù)值解關(guān)于時(shí)間步長(zhǎng)h的收斂行為。

表1 算例1中α=0.1，0.4，最大誤差隨時(shí)間步長(zhǎng)h的變化與收斂階Tab.1 Maximum errors and decay rate as functions of h with α=0.1，0.4，for Example 1

表2 算例1中α=0.7，0.99，最大誤差隨時(shí)間步長(zhǎng)h的變化與收斂階Tab.2 Maximum errors and decay rate as functions of h withα=0.7，0.99，for Example 1

例2，考慮初值問(wèn)題(1)-(2):

其中λ是一個(gè)常數(shù)。注意到此時(shí)f是y的線性函數(shù)，相應(yīng)的精確解為y(x)=x3+α。

首先，取λ=1，表3、表4重復(fù)例1的計(jì)算過(guò)程，我們得到取α從0.1到0.99時(shí)，最大誤差接近于3-α，這支持了理論分析的結(jié)果。

表4 算例2中α=0.7，0.99，最大誤差隨時(shí)間步長(zhǎng)h的變化與收斂階Tab.4 Maximum errors and decay rate as functions of h withα=0.7，0.99 for Example 2