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

(南京航空航天大學 機械結(jié)構(gòu)力學及控制國家重點實驗室，南京 210016)

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

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

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

1 分數(shù)階導數(shù)的記憶性

Grünwald-Letnikov定義為

(1)

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

(2)

(3)

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

(4)

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

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

2 改進的短記憶原理

Grünwald-Letnikov定義下分數(shù)階導數(shù)計算時取值如圖1所示。計算第k個時間步時，需要考慮前面所有時刻的影響。如果需要計算的時間步很多時，計算量將是非常龐大。此外為了獲得較為精確的數(shù)值解，通常以減小時間步長為手段，這是計算量增大另一個因素。因此有了短記憶原理的提出，如圖2所示。L為截取的記憶時間，反映在公式上是截斷Nt項。如果步長為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ù)項數(shù)的截斷。當計算時間超過L時，則把超出的部分時刻取值步長調(diào)整為初始步長的倍數(shù)。調(diào)整步長后，用小步長計算的時刻不用再次計算，只計算小步長范圍外的部分，如圖3、圖4和圖5所示。圖中二項式截斷項數(shù)為Nt=10，步長增大倍數(shù)為n=2。

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

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

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

圖3中計算步數(shù)k≤Nt，未到截斷條件，則按照原始定義計算。表達式為

(6)

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

(7)

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

(8)

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

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

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

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

3 算例分析

3.1 分數(shù)階阻尼受迫振動

采用帶有分數(shù)階阻尼的振動微分方程為例，方程形式為

(9)

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

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

圖8 改進的短記憶原理方法(步長同為h=0.000 1，放大倍數(shù)同為n=10，不同的截斷項數(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 改進的短記憶原理方法(步長同為h=0.000 1，截斷項數(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 分數(shù)階非線性Duffing方程

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

(10)

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

圖10 分數(shù)階Duffing方程的位移圖和相圖分別用龍格庫塔法和改進的短記憶原理方法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 分數(shù)階Lorenz混沌系統(tǒng)

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

(11)

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

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

4 結(jié) 論

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