馬瑞群， 張 波， 員海瑋， 韓景龍

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

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

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

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

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

Grünwald-Letnikov定義為

(1)

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

(2)

(3)

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

(4)

wj是從當(dāng)前時(shí)刻往前數(shù)j個(gè)步長(zhǎng)間隔的時(shí)刻的函數(shù)值的加權(quán)系數(shù)，也即二項(xiàng)式系數(shù)。從遞推公式的等號(hào)左邊的系數(shù)小于1可知這個(gè)加權(quán)系數(shù)隨著時(shí)間間隔的增大而減小。表1給出了當(dāng)α=0.5時(shí)的部分二項(xiàng)式系數(shù)?？梢钥吹郊訖?quán)系數(shù)的絕對(duì)值在開(kāi)始階段迅速減小，隨后減小的速度越來(lái)越慢。從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萬(wàn)項(xiàng)，雖然每一個(gè)時(shí)刻對(duì)當(dāng)前時(shí)刻影響很小，但是項(xiàng)數(shù)如此之多依然會(huì)帶來(lái)很大的誤差。

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

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

Grünwald-Letnikov定義下分?jǐn)?shù)階導(dǎo)數(shù)計(jì)算時(shí)取值如圖1所示。計(jì)算第k個(gè)時(shí)間步時(shí)，需要考慮前面所有時(shí)刻的影響。如果需要計(jì)算的時(shí)間步很多時(shí)，計(jì)算量將是非常龐大。此外為了獲得較為精確的數(shù)值解，通常以減小時(shí)間步長(zhǎng)為手段，這是計(jì)算量增大另一個(gè)因素。因此有了短記憶原理的提出，如圖2所示。L為截取的記憶時(shí)間，反映在公式上是截?cái)郚t項(xiàng)。如果步長(zhǎ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ù)是，隨著時(shí)間的推移，久遠(yuǎn)的時(shí)刻對(duì)當(dāng)前時(shí)刻的影響越來(lái)越小，因此就進(jìn)行有效截?cái)嗵幚?。但是這種截?cái)嗤鶗?huì)帶來(lái)不可忽視的誤差。所以對(duì)短記憶原理進(jìn)行改進(jìn)。改進(jìn)的思路是把短記憶原理對(duì)時(shí)間的截?cái)嗵鎿Q為對(duì)二項(xiàng)式系數(shù)項(xiàng)數(shù)的截?cái)?。?dāng)計(jì)算時(shí)間超過(guò)L時(shí)，則把超出的部分時(shí)刻取值步長(zhǎng)調(diào)整為初始步長(zhǎng)的倍數(shù)。調(diào)整步長(zhǎng)后，用小步長(zhǎng)計(jì)算的時(shí)刻不用再次計(jì)算，只計(jì)算小步長(zhǎng)范圍外的部分，如圖3、圖4和圖5所示。圖中二項(xiàng)式截?cái)囗?xiàng)數(shù)為Nt=10，步長(zhǎng)增大倍數(shù)為n=2。

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

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

圖5 計(jì)算時(shí)間步數(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]。以此類(lèi)推，可以得到更多的步長(zhǎng)表達(dá)式。

理論上，初始截?cái)囗?xiàng)數(shù)Nt越大，所取得的結(jié)果越精確。但是考慮計(jì)算成本，會(huì)選擇一個(gè)比較合適的值。至于步長(zhǎng)放大倍數(shù)，換一個(gè)角度考慮。步長(zhǎng)放大后，所取函數(shù)值可以看作為該間隔的代表。式(7)等號(hào)右邊可以把公因子1/h1α提出來(lái)，則第二項(xiàng)包含二項(xiàng)式系數(shù)的第i項(xiàng)影響系數(shù)為(1/nα)wi，其與未放大步長(zhǎng)時(shí)該函數(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í)的部分影響系數(shù)Tab.2 Partial influence coefficient at α=0.5

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

3 算例分析

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

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

(9)

接下來(lái)考慮不同截?cái)囗?xiàng)數(shù)的情況。設(shè)置所有的基礎(chǔ)步長(zhǎng)h=0.000 1，所有步長(zhǎng)放大倍數(shù)為n=10。斷項(xiàng)數(shù)分別為Nt=3 000，Nt=10 000，Nt=30 000，所對(duì)應(yīng)的初始記憶時(shí)間為0.3 s、1 s和3 s。具體數(shù)值結(jié)果與圖6無(wú)異，由于誤差量級(jí)較小，看不到區(qū)別。其誤差結(jié)果如圖8所示。黑色實(shí)線為Nt=3 000時(shí)與原始定義數(shù)值解之間的誤差，與另外兩條線相比幅值比較大，是因?yàn)槌跏加洃洉r(shí)間較短造成的。虛線Nt=10 000與點(diǎn)劃線Nt=30 000能很好地逼近原始定義數(shù)值解。這說(shuō)明初始記憶時(shí)間越長(zhǎng)所取得的數(shù)值解越精確。從計(jì)算量上來(lái)看，Nt=3 000，Nt=10 000，Nt=30 000時(shí)，最后一步參與計(jì)算的時(shí)刻點(diǎn)數(shù)分別為8 600項(xiàng)、23 000項(xiàng)和59 000項(xiàng)?？梢钥吹饺N情況計(jì)算量均比原始定義方法少得多，另外也可以得到初始記憶時(shí)間越長(zhǎng)所得結(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給出了相同的步長(zhǎng)和截?cái)囗?xiàng)數(shù)，不同的的步長(zhǎng)放大倍數(shù)下，數(shù)值解與原始定義之間的誤差情況。其中實(shí)線、虛線和點(diǎn)劃線分別為放大5倍、10倍和20倍的結(jié)果與原始定義數(shù)值解的誤差。可以看到在放大20倍時(shí)誤差較大，這與表二的結(jié)果吻合。另外，放大5倍與放大10倍時(shí)相比，誤差并沒(méi)有相應(yīng)的減少太多。這從表二中也能看出來(lái)，兩次放大的影響系數(shù)與各自對(duì)應(yīng)的放大倍數(shù)都很接近。計(jì)算量上，n=5，n=10，n=20時(shí)最后一個(gè)時(shí)間步參與計(jì)算的時(shí)刻點(diǎn)數(shù)分別為10 400項(xiàng)、8 600項(xiàng)以及6 950項(xiàng)。放大5倍的計(jì)算量比放大10倍大了不少，但是精度卻并沒(méi)有提升太多。但是放大10倍與放大20倍相比，其所增加的計(jì)算量和誤差的減小是符合我們期望的。

圖8 改進(jìn)的短記憶原理方法(步長(zhǎng)同為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)的短記憶原理方法(步長(zhǎng)同為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í)動(dòng)力系統(tǒng)幾乎總是含有各種各樣的非線性因素。Duffing方程是機(jī)械振動(dòng)學(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，時(shí)間步長(zhǎng)均取h=0.001，計(jì)算1 000 s，截取后700 s。式(10)分別采用4階龍格庫(kù)塔法和改進(jìn)的短記憶原理方法進(jìn)行求解。龍格庫(kù)塔法的位移曲線如圖10(a)所示。短記憶原理方法的截?cái)囗?xiàng)數(shù)為Nt=10 000，步長(zhǎng)放大倍數(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方程的位移圖和相圖分別用龍格庫(kù)塔法和改進(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)是一個(gè)典型的多自由度混沌系統(tǒng)。簡(jiǎn)化的Lorenz系統(tǒng)只含一個(gè)系統(tǒng)參數(shù)，其方程為

(11)

系統(tǒng)參數(shù)設(shè)定為c=10。分?jǐn)?shù)階數(shù)α，β，γ均取0.98。用預(yù)估校正法與改進(jìn)的短記憶原理方法計(jì)算50 s，時(shí)間步長(zhǎng)均取h=0.001，取后40 s數(shù)值結(jié)果進(jìn)行比較，如圖11所示?？梢钥吹?，改進(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ì)算中通常采取減小步長(zhǎng)來(lái)獲取高精度數(shù)值解，而分?jǐn)?shù)階導(dǎo)數(shù)具有記憶性，每一個(gè)時(shí)刻都與前面所有時(shí)刻相關(guān)，所以減小步長(zhǎng)的同時(shí)會(huì)帶來(lái)巨大的計(jì)算量。經(jīng)典的短記憶原理，在減少計(jì)算量上有很大貢獻(xiàn)，但是造成的誤差不可忽視。本文提出一種基于Grünwald-Letnikov定義的改進(jìn)的短記憶原理方法，并應(yīng)用于分?jǐn)?shù)階阻尼振動(dòng)方程、分?jǐn)?shù)階非線性Duffing方程和分?jǐn)?shù)階Lorenz混沌系統(tǒng)。分?jǐn)?shù)階阻尼振動(dòng)方程的解通過(guò)與經(jīng)典的短記憶原理方法相比，改進(jìn)的短記憶原理方法能極大的減小誤差，并且沒(méi)有增加計(jì)算量。通過(guò)分?jǐn)?shù)階非線性Duffing方程的算例，說(shuō)明該方法也適用于模擬非線性系統(tǒng)。改進(jìn)的短記憶原理方法應(yīng)用于Lorenz混沌系統(tǒng)也有很好的數(shù)值結(jié)果。在改進(jìn)的短記憶原理方法中，即使減小步長(zhǎng)提高精度，也不會(huì)帶來(lái)很大的計(jì)算量的變化。在使用改進(jìn)的短記憶原理方法時(shí)，截?cái)囗?xiàng)數(shù)在計(jì)算機(jī)和時(shí)間允許的情況下盡可能取得更大，而放大倍數(shù)的選取可通過(guò)計(jì)算得到的影響系數(shù)與放大倍數(shù)匹配。