馬瑞群， 張 波， 員?，|， 韓景龍

(南京航空航天大學(xué) 機(jī)械結(jié)構(gòu)力學(xué)及控制國家重點(diǎn)實(shí)驗(yàn)室，南京 210016)

分?jǐn)?shù)階導(dǎo)數(shù)的概念可以追溯到微分學(xué)開始的時候，1695年，Leibniz[1]和L′Hpital在往來書信中開始討論分?jǐn)?shù)階導(dǎo)數(shù)。隨后，通過Euler、Liouville、Riemann、Letnikov和Coputo等偉大數(shù)學(xué)家的共同努力，奠定了分?jǐn)?shù)導(dǎo)數(shù)的數(shù)學(xué)基礎(chǔ)。然而，直到最近，分?jǐn)?shù)階導(dǎo)數(shù)在科學(xué)和工程的各個分支中的重要應(yīng)用才得以確立。分?jǐn)?shù)階微分方程在物理學(xué)[2-5]、化學(xué)[6-7]和工程學(xué)[8-10]中扮演著越來越重要的角色。從時間上而言，整數(shù)階導(dǎo)數(shù)所表征的是物理過程某時刻的物理量的變化，而分?jǐn)?shù)階導(dǎo)數(shù)所表征的性質(zhì)則與該現(xiàn)象的整個發(fā)展歷史有關(guān)。整數(shù)階空間導(dǎo)數(shù)描述的是一個物理過程在空間中某一確定位置的局部性質(zhì)，而分?jǐn)?shù)階導(dǎo)數(shù)所描述的性質(zhì)則與物理過程所涉及的整個空間有關(guān)。盡管分?jǐn)?shù)階導(dǎo)數(shù)有很多種定義，但最常用的定義是Grünwald-Letnikov、Riemann-Liouville和Caputo分?jǐn)?shù)導(dǎo)數(shù)。

由于時間分?jǐn)?shù)階導(dǎo)數(shù)具有的時間記憶性，其在數(shù)值求解過程中為了得到更加精確的結(jié)果會有大量的數(shù)據(jù)參與計(jì)算。分?jǐn)?shù)階微分方程數(shù)值近似解的方法已經(jīng)被廣泛研究[11-17]，但是通常在計(jì)算效率、復(fù)雜性和結(jié)果近似的精度之間總是存在權(quán)衡。文獻(xiàn)[18]提出了一種沒有記憶效應(yīng)的高效仿真方法。文獻(xiàn)[19]提出短記憶原理(short memory principle，SMP)，是指截取最近的一個時間段，而對影響比較小且較遠(yuǎn)的時刻選擇忽略。文獻(xiàn)[20]利用短記憶原理驗(yàn)證了周期函數(shù)的分?jǐn)?shù)階導(dǎo)數(shù)依然是周期函數(shù)。文獻(xiàn)[21]結(jié)合短記憶原理和Grünwald-Letnikov定義，研究了分?jǐn)?shù)階系統(tǒng)的模型參考自適應(yīng)控制。Wei等[22]根據(jù)Grünwald-Letnikov定義下的經(jīng)典短記憶原理，提出并研究了幾種新穎的短時記憶原理。傳統(tǒng)的短記憶原理可能會帶來較大的誤差，例如自由振動不能回到平衡位置。尤其當(dāng)分?jǐn)?shù)階數(shù)α趨于0時，為了保證計(jì)算精度用短記憶原理截取的記憶時間可能極其巨大。文獻(xiàn)[23-25]提到用嵌套網(wǎng)格方法解決長時記憶，其中：Ford等提出在時間段的兩端用小步長，但是隨著時間推進(jìn)，初始時刻附近的系數(shù)將變得很小且每一步都需要重新計(jì)算；Diethelm等所提方法則同樣需要每一步重新計(jì)算系數(shù)。

本文提出一種基于Grünwald-Letnikov定義改進(jìn)的短記憶原理。這種改進(jìn)是把對記憶時間的截?cái)嗾{(diào)整為對二項(xiàng)式系數(shù)的截?cái)?。為了保證計(jì)算的精度，初始采用小步長。如果計(jì)算過的數(shù)據(jù)點(diǎn)數(shù)超過二項(xiàng)式系數(shù)的項(xiàng)數(shù)，則增大步長使有限的二項(xiàng)式系數(shù)能覆蓋更大的時間區(qū)域。這種改進(jìn)的目的是用盡量少的數(shù)據(jù)量求得比較精確的結(jié)果。文章中用單自由度分?jǐn)?shù)階阻尼受迫振動算例驗(yàn)證了方法的準(zhǔn)確性和可靠性。

1 分?jǐn)?shù)階導(dǎo)數(shù)的記憶性

Grünwald-Letnikov定義為

(1)

式中，h為計(jì)算步長。如果選擇的計(jì)算步長足夠小，則求極限操作可以忽略，變?yōu)?/p>

(2)

(3)

為了避免Gamma函數(shù)的計(jì)算，采用遞推形式

(4)

wj是從當(dāng)前時刻往前數(shù)j個步長間隔的時刻的函數(shù)值的加權(quán)系數(shù)，也即二項(xiàng)式系數(shù)。從遞推公式的等號左邊的系數(shù)小于1可知這個加權(quán)系數(shù)隨著時間間隔的增大而減小。表1給出了當(dāng)α=0.5時的部分二項(xiàng)式系數(shù)?？梢钥吹郊訖?quán)系數(shù)的絕對值在開始階段迅速減小，隨后減小的速度越來越慢。從1～-0.01只計(jì)算了9項(xiàng)，后面依次減小10倍的項(xiàng)分別為w43，w200，w926，w4301，w19965，w92668和w430126。從-1.00×10-8～-1.00×10-9需要經(jīng)歷將近34萬項(xiàng)，雖然每一個時刻對當(dāng)前時刻影響很小，但是項(xiàng)數(shù)如此之多依然會帶來很大的誤差。

表1 當(dāng)α=0.5時的部分二項(xiàng)式系數(shù)Tab.1 Partial binomial coefficient at α=0.5

2 改進(jìn)的短記憶原理

Grünwald-Letnikov定義下分?jǐn)?shù)階導(dǎo)數(shù)計(jì)算時取值如圖1所示。計(jì)算第k個時間步時，需要考慮前面所有時刻的影響。如果需要計(jì)算的時間步很多時，計(jì)算量將是非常龐大。此外為了獲得較為精確的數(shù)值解，通常以減小時間步長為手段，這是計(jì)算量增大另一個因素。因此有了短記憶原理的提出，如圖2所示。L為截取的記憶時間，反映在公式上是截?cái)郚t項(xiàng)。如果步長為h，則L=Nt×h。短記憶原理公式為

(5)

圖1 Grünwald-Letnikov原始定義示意圖Fig.1 Grünwald-Letnikov original definition diagram

圖2 短記憶原理示意圖Fig.2 Schematic diagram of short memory

短記憶原理的提出的根據(jù)是，隨著時間的推移，久遠(yuǎn)的時刻對當(dāng)前時刻的影響越來越小，因此就進(jìn)行有效截?cái)嗵幚?。但是這種截?cái)嗤鶗聿豢珊鲆暤恼`差。所以對短記憶原理進(jìn)行改進(jìn)。改進(jìn)的思路是把短記憶原理對時間的截?cái)嗵鎿Q為對二項(xiàng)式系數(shù)項(xiàng)數(shù)的截?cái)?。?dāng)計(jì)算時間超過L時，則把超出的部分時刻取值步長調(diào)整為初始步長的倍數(shù)。調(diào)整步長后，用小步長計(jì)算的時刻不用再次計(jì)算，只計(jì)算小步長范圍外的部分，如圖3、圖4和圖5所示。圖中二項(xiàng)式截?cái)囗?xiàng)數(shù)為Nt=10，步長增大倍數(shù)為n=2。

圖3 計(jì)算時間步數(shù)k≤NtFig.3 Calculation time steps k≤Nt

圖4 計(jì)算時間步數(shù)Nt

圖5 計(jì)算時間步數(shù)nNtn2NtFig.5 Calculate the number of time steps nNtn2Nt

圖3中計(jì)算步數(shù)k≤Nt，未到截?cái)鄺l件，則按照原始定義計(jì)算。表達(dá)式為

(6)

圖4中計(jì)算步數(shù)Nt

(7)

式中：a,b分別為[Nt/n]+1，[k/n];[·]為取整。步數(shù)nNt

(8)

式中，d為[k/n2]。以此類推，可以得到更多的步長表達(dá)式。

理論上，初始截?cái)囗?xiàng)數(shù)Nt越大，所取得的結(jié)果越精確。但是考慮計(jì)算成本，會選擇一個比較合適的值。至于步長放大倍數(shù)，換一個角度考慮。步長放大后，所取函數(shù)值可以看作為該間隔的代表。式(7)等號右邊可以把公因子1/h1α提出來，則第二項(xiàng)包含二項(xiàng)式系數(shù)的第i項(xiàng)影響系數(shù)為(1/nα)wi，其與未放大步長時該函數(shù)值的影響系數(shù)比值幾乎等于放大倍數(shù)n。以截?cái)囗?xiàng)數(shù)Nt=100，放大倍數(shù)分別為n=5，n=10，n=20為例，α=0.5部分比值如表2所示。

表2 當(dāng)α=0.5時的部分影響系數(shù)Tab.2 Partial influence coefficient at α=0.5

可以看到第一次放大后，影響系數(shù)隨著步數(shù)遞增逐漸減小，越來越接近放大倍數(shù)n。隨著放大倍數(shù)的增大，第一個影響系數(shù)變大，尤其當(dāng)n=20時，大過該步長內(nèi)初始步長的個數(shù)，因此會帶來較大誤差。另外所取時刻的函數(shù)值在該間隔內(nèi)是否具有代表性，即所取函數(shù)值是否更接近該間隔內(nèi)函數(shù)值的平均值。因?yàn)閿U(kuò)大步長取函數(shù)值是一個動態(tài)的遍歷過程，所以在長時計(jì)算中引起的誤差也有限。綜上所述，在計(jì)算條件允許下盡量截取更多的二項(xiàng)式系數(shù)，以及選取較小的步長放大倍數(shù)，才能獲得更為精確的數(shù)值解。

3 算例分析

3.1 分?jǐn)?shù)階阻尼受迫振動

采用帶有分?jǐn)?shù)階阻尼的振動微分方程為例，方程形式為

(9)

接下來考慮不同截?cái)囗?xiàng)數(shù)的情況。設(shè)置所有的基礎(chǔ)步長h=0.000 1，所有步長放大倍數(shù)為n=10。斷項(xiàng)數(shù)分別為Nt=3 000，Nt=10 000，Nt=30 000，所對應(yīng)的初始記憶時間為0.3 s、1 s和3 s。具體數(shù)值結(jié)果與圖6無異，由于誤差量級較小，看不到區(qū)別。其誤差結(jié)果如圖8所示。黑色實(shí)線為Nt=3 000時與原始定義數(shù)值解之間的誤差，與另外兩條線相比幅值比較大，是因?yàn)槌跏加洃洉r間較短造成的。虛線Nt=10 000與點(diǎn)劃線Nt=30 000能很好地逼近原始定義數(shù)值解。這說明初始記憶時間越長所取得的數(shù)值解越精確。從計(jì)算量上來看，Nt=3 000，Nt=10 000，Nt=30 000時，最后一步參與計(jì)算的時刻點(diǎn)數(shù)分別為8 600項(xiàng)、23 000項(xiàng)和59 000項(xiàng)。可以看到三種情況計(jì)算量均比原始定義方法少得多，另外也可以得到初始記憶時間越長所得結(jié)果越精確的結(jié)論。

圖6 三種數(shù)值解(原始定義、短記憶原理和改進(jìn)的短記憶原理)Fig.6 Three numerical solutions (original definition, short memory principle and improved short memory principle)

圖7 兩種短記憶原理方法與原始定義方法間的誤差Fig.7 The error between the two short memory principle methods and the original definition method

圖9給出了相同的步長和截?cái)囗?xiàng)數(shù)，不同的的步長放大倍數(shù)下，數(shù)值解與原始定義之間的誤差情況。其中實(shí)線、虛線和點(diǎn)劃線分別為放大5倍、10倍和20倍的結(jié)果與原始定義數(shù)值解的誤差?？梢钥吹皆诜糯?0倍時誤差較大，這與表二的結(jié)果吻合。另外，放大5倍與放大10倍時相比，誤差并沒有相應(yīng)的減少太多。這從表二中也能看出來，兩次放大的影響系數(shù)與各自對應(yīng)的放大倍數(shù)都很接近。計(jì)算量上，n=5，n=10，n=20時最后一個時間步參與計(jì)算的時刻點(diǎn)數(shù)分別為10 400項(xiàng)、8 600項(xiàng)以及6 950項(xiàng)。放大5倍的計(jì)算量比放大10倍大了不少，但是精度卻并沒有提升太多。但是放大10倍與放大20倍相比，其所增加的計(jì)算量和誤差的減小是符合我們期望的。

圖8 改進(jìn)的短記憶原理方法(步長同為h=0.000 1，放大倍數(shù)同為n=10，不同的截?cái)囗?xiàng)數(shù)Nt)與原始定義方法間數(shù)值解的誤差Fig.8 The error of the numerical solution between the improvedshort memory principle method (the step size h=0.000 1 and magnification factor n=10 are the same, and the number of truncation items Nt are different) and the original definition method

圖9 改進(jìn)的短記憶原理方法(步長同為h=0.000 1，截?cái)囗?xiàng)數(shù)Nt=3 000，不同的放大倍數(shù)n)與原始定義方法間數(shù)值解的誤差Fig.9 The error between the improved short memory principle method (the step size h=0.000 1 and the number of truncation items Nt=3 000 are the same，and magnification factors n are different) and the original definition method

3.2 分?jǐn)?shù)階非線性Duffing方程

工程中的真實(shí)動力系統(tǒng)幾乎總是含有各種各樣的非線性因素。Duffing方程是機(jī)械振動學(xué)、物理學(xué)、生物學(xué)和神經(jīng)學(xué)等領(lǐng)域廣泛應(yīng)用的數(shù)學(xué)模型。具有分?jǐn)?shù)階阻尼的Duffing系統(tǒng)為

(10)

分?jǐn)?shù)階數(shù)α=0.5，時間步長均取h=0.001，計(jì)算1 000 s，截取后700 s。式(10)分別采用4階龍格庫塔法和改進(jìn)的短記憶原理方法進(jìn)行求解。龍格庫塔法的位移曲線如圖10(a)所示。短記憶原理方法的截?cái)囗?xiàng)數(shù)為Nt=10 000，步長放大倍數(shù)為n=5。計(jì)算結(jié)果的位移曲線如圖10(b)所示。從位移圖上看兩者的結(jié)果很相似，然后再看兩種方法的系統(tǒng)相圖，如圖10(c)和圖10(d)所示。從位移圖和相圖上看，該改進(jìn)的短記憶原理方法能很好的模擬分?jǐn)?shù)階非線性系統(tǒng)。

圖10 分?jǐn)?shù)階Duffing方程的位移圖和相圖分別用龍格庫塔法和改進(jìn)的短記憶原理方法Fig.10 The displacement diagram and phase diagram of the fractional Duffing equation use Runge-Kutta method and improved short memory principle method respectively

3.3 分?jǐn)?shù)階Lorenz混沌系統(tǒng)

分?jǐn)?shù)階Lorenz系統(tǒng)是一個典型的多自由度混沌系統(tǒng)。簡化的Lorenz系統(tǒng)只含一個系統(tǒng)參數(shù)，其方程為

(11)

系統(tǒng)參數(shù)設(shè)定為c=10。分?jǐn)?shù)階數(shù)α，β，γ均取0.98。用預(yù)估校正法與改進(jìn)的短記憶原理方法計(jì)算50 s，時間步長均取h=0.001，取后40 s數(shù)值結(jié)果進(jìn)行比較，如圖11所示?？梢钥吹剑倪M(jìn)的短記憶原理方法能很好地模擬多自由度系統(tǒng)。

圖11 分?jǐn)?shù)階Lorenz混沌吸引子部分相圖預(yù)估校正法和改進(jìn)的短記憶原理方法Fig.11 Fractional Lorenz chaotic attractor partial phase diagram prediction correction method and improved short memory principle method

4 結(jié) 論

在數(shù)值計(jì)算中通常采取減小步長來獲取高精度數(shù)值解，而分?jǐn)?shù)階導(dǎo)數(shù)具有記憶性，每一個時刻都與前面所有時刻相關(guān)，所以減小步長的同時會帶來巨大的計(jì)算量。經(jīng)典的短記憶原理，在減少計(jì)算量上有很大貢獻(xiàn)，但是造成的誤差不可忽視。本文提出一種基于Grünwald-Letnikov定義的改進(jìn)的短記憶原理方法，并應(yīng)用于分?jǐn)?shù)階阻尼振動方程、分?jǐn)?shù)階非線性Duffing方程和分?jǐn)?shù)階Lorenz混沌系統(tǒng)。分?jǐn)?shù)階阻尼振動方程的解通過與經(jīng)典的短記憶原理方法相比，改進(jìn)的短記憶原理方法能極大的減小誤差，并且沒有增加計(jì)算量。通過分?jǐn)?shù)階非線性Duffing方程的算例，說明該方法也適用于模擬非線性系統(tǒng)。改進(jìn)的短記憶原理方法應(yīng)用于Lorenz混沌系統(tǒng)也有很好的數(shù)值結(jié)果。在改進(jìn)的短記憶原理方法中，即使減小步長提高精度，也不會帶來很大的計(jì)算量的變化。在使用改進(jìn)的短記憶原理方法時，截?cái)囗?xiàng)數(shù)在計(jì)算機(jī)和時間允許的情況下盡可能取得更大，而放大倍數(shù)的選取可通過計(jì)算得到的影響系數(shù)與放大倍數(shù)匹配。