姚偉岸, 高 強, 譚述君, 吳 鋒

(大連理工大學(xué) 力學(xué)與航空航天學(xué)院 工業(yè)裝備結(jié)構(gòu)分析優(yōu)化與CAE軟件全國重點實驗室, 大連 116024)

1 引 言

精細積分方法的基本思想由鐘萬勰院士1991年首次提出[1],并很快應(yīng)用到結(jié)構(gòu)動力初值問題的求解[2]。隨后,結(jié)合計算結(jié)構(gòu)力學(xué)與最優(yōu)控制的模擬理論,鐘萬勰院士將精細積分方法的思想推廣到兩點邊值問題[3],并持續(xù)不斷推進精細積分方法的研究工作[4]。精細積分方法提出之后,經(jīng)過鐘萬勰院士眾多弟子及不同領(lǐng)域?qū)W術(shù)同行的不斷發(fā)展,在非齊次方程精細模擬、時變動力問題、非線性動力問題、復(fù)雜大規(guī)模問題和特殊矩陣函數(shù)等基本理論和算法方面不斷取得研究成果,并已經(jīng)拓展應(yīng)用到隨機動力響應(yīng)分析、熱傳導(dǎo)分析、動力彈塑性硬化和軟化問題、動態(tài)載荷識別、周期結(jié)構(gòu)能帶分析、最優(yōu)控制問題、Floquet轉(zhuǎn)移矩陣計算、對稱和非對稱Riccati微分方程求解和Maxwell方程求解等眾多領(lǐng)域。

精細積分方法是一種區(qū)別于傳統(tǒng)差分類數(shù)值方法的全新的初值問題數(shù)值積分方法,其突出特色體現(xiàn)在:(1) 數(shù)值積分過程中,基于加法定理和2N類算法,通過僅計算和存儲迭代過程中的增量,將計算過程中的截斷和舍入誤差降低到計算機精度之外,從而能夠給出計算機浮點意義上的精確解,并能在合理積分步長范圍內(nèi)不發(fā)生穩(wěn)定性和剛性問題; (2) 精細積分方法不僅是顯式格式,并且是無條件穩(wěn)定的,突破了差分類算法只有隱式格式才具有無條件穩(wěn)定性的限制; (3) 精細積分方法具有零振幅衰減率、零周期擴散率和無超越性等優(yōu)點。由于算法的優(yōu)異性能,能夠獲得遠優(yōu)于傳統(tǒng)差分算法的數(shù)值結(jié)果,特別是在處理剛性問題時具有差分類算法無法比擬的穩(wěn)定性和精度。

本文作者在跟隨鐘萬勰院士讀博士和后期工作期間,有幸在鐘先生的教誨和指導(dǎo)下從不同方面參與到精細積分方法的理論研究和應(yīng)用工作,并不斷跟蹤精細積分方法的發(fā)展情況[5]。本文對精細積分方法的基本思想、深入發(fā)展和在各領(lǐng)域的應(yīng)用做一系統(tǒng)綜述。

2 精細積分方法的基本思想

考慮常微分方程矩陣-向量形式的一般性表達

(1)

當系統(tǒng)矩陣A和向量f與狀態(tài)v不相關(guān)時,按線性微分方程的求解理論[6],狀態(tài)微分方程(1)的通解可以由狀態(tài)傳遞矩陣和Duhamel積分表達為

(2)

式中Φ(t,τ)為狀態(tài)傳遞矩陣,且具有以下性質(zhì)

(1)Φ(t,t)=I

(3)

(2)Φ(t,t1)=Φ(t,t2)Φ(t2,t1)

(4)

(5)

式中 用A(t)描述表明這對于時變系統(tǒng)也適用。

式(2)表明,狀態(tài)傳遞矩陣Φ是求解線性微分方程的關(guān)鍵。對于時不變系統(tǒng),A為常數(shù)矩陣時,狀態(tài)傳遞矩陣僅與區(qū)段長度有關(guān),即

Φ(t,t1)=Φ(t-t1)=eA(t-t1)

(6)

式中 矩陣指數(shù)

(7)

雖然從數(shù)學(xué)上,矩陣指數(shù)可采用本征向量展開方法給出,但在實際數(shù)值計算中并不是十分有效,尤其出現(xiàn)或接近出現(xiàn)Jordan型本征解時。關(guān)于矩陣指數(shù)已經(jīng)提出了很多計算方法,但仍不夠理想。文獻[7]給出了19種可疑的算法,且25年后再予以回顧[8],指出問題并未完全解決。

2.1 精細積分方法的提出

鐘萬勰院士研究連續(xù)時間LQ控制本征對求解時,針對其核心問題即矩陣指數(shù)的計算,首次提出了精細積分方法[1],并達到計算機意義上的精確解。本小節(jié)綜合文獻[1,2,9]的論述,介紹精細積分方法的基本思想。

(8)

根據(jù)矩陣指數(shù)的加法定理可給出

(9)

式中m=2N,而N為任意正整數(shù),如可選N=20,則m=1048576。亦稱為2N類算法。

eAτ≈I+Aτ+(Aτ)2/2+…+(Aτ)q/q!=I+Ta

(10)

式中Ta陣是一個小量的矩陣,為

Ta=Aτ+(Aτ)2/2+…+(Aτ)q/q!

(11)

計算中至關(guān)重要的一點是矩陣指數(shù)的存儲只能是增量Ta,而不是(I+Ta)。因為Ta很小,當其與矩陣I相加時,就會成為其尾數(shù),在計算機的舍入操作中,其精度將喪失殆盡。為了提高計算T陣的效率,運用2N類算法先將式(9)作分解

(12)

這種分解一直做下去,共N次。

其次應(yīng)注意,對任意矩陣Tb和Tc有

(I+Tb)(I+Tc)≡I+Tb+Tc+TbTc

(13)

當Tb和Tc很小時,不應(yīng)加上I之后再執(zhí)行乘法。因此式(12)的N次乘法相當于

(14)

T=I+Ta

(15)

由于N次乘法后,Ta已不再是很小的矩陣,這個加法已沒有嚴重的舍入誤差。

以上是矩陣指數(shù)計算的精細積分方法的思想。文獻[1]指出,精細積分有兩個要點,(1) 運用矩陣指數(shù)的加法定理,即2N類算法; (2) 注意力放在增量上,而不是其全量。

2.2 基于Padé級數(shù)近似的精細積分方法

Rpq(Aτ)=[Dpq(Aτ)]-1Npq(Aτ)

(16)

式中

(17)

(18)

顯然,當q=0時,Rp0(Aτ)是eAτ的p階Taylor近似,即Taylor近似是Padé近似的特例。在實際應(yīng)用中通常令p=q,即q階對角Padé近似。下面介紹的Padé近似均表示對角Padé近似。

參照精細積分方法,將精細區(qū)段矩陣指數(shù)的q階Padé近似的增量分離出來表示為

Rqq(Aτ)=I+Ra

(19)

同樣Nqq(Aτ)和Dqq(Aτ)的增量也分離出來

Nqq(Aτ)=I+Na,Dqq(Aτ)=I+Da

(20)

式中

(21)

將式(19,20)代入式(16),可以導(dǎo)出

Ra=(I+Da)-1(Na-Da)

(22)

獲得了Padé近似的增量部分Ra后,再采用與2.1節(jié)相同的方法,結(jié)合矩陣指數(shù)加法定理,并在執(zhí)行過程中只存儲和計算增量部分,就得到了基于Padé級數(shù)近似的矩陣指數(shù)精細積分方法。

由于在精細區(qū)段τ上采用了精度更高的Padé級數(shù)近似,因此在選擇相同參數(shù)(N,q)的情況下,基于Padé級數(shù)近似的精細積分方法比基于Taylor級數(shù)近似的精細積分方法計算精度更高和穩(wěn)定性更好。但是Padé級數(shù)增量部分Ra涉及矩陣求逆,這一定程度上影響了其計算效率,尤其是矩陣維數(shù)較大時。進一步需要指出的是,如矩陣A是Hamilton矩陣,則其對應(yīng)的矩陣指數(shù)是辛矩陣。此情況可以證明基于Padé近似的精細積分方法給出的矩陣指數(shù)是辛矩陣,因此是保辛算法[11,12]。

2.3 誤差分析與算法設(shè)計

無論是基于Taylor級數(shù)近似還是基于Padé級數(shù)近似的精細積分方法,其精細區(qū)段劃分系數(shù)N和截斷項數(shù)q的選擇直接影響到矩陣指數(shù)數(shù)值逼近的精度,也直接影響到計算效率。文獻[13,14]分析了基于Taylor級數(shù)逼近的精細積分方法的誤差估計,并給出了(N,q)選擇的經(jīng)驗及一些建議。文獻[10]研究了當展開項數(shù)q=4時,系數(shù)N的選擇,而文獻[15]則分析了Hamilton矩陣指數(shù)采用一階和二階甚至任意階Padé級數(shù)近似的誤差估計等。上面對(N,q)的選擇是在給定級數(shù)截斷項數(shù)q(主要憑經(jīng)驗)的情況下,選擇系數(shù)N以保證其計算精度,但這樣得到的(N,q)組合不能保證其計算量最小。

文獻[16]進一步基于誤差分析理論,推導(dǎo)出給定精度下采用Padé級數(shù)逼近時參數(shù)(N,q)滿足的條件?？紤]到矩陣指數(shù)精細積分方法的計算量主要為(N+q-1)次矩陣乘法,因此可將參數(shù)(N,q)的選擇描述為下面優(yōu)化問題的求解

(23)

式中 EPS為指定精度

(24)

進一步,文獻[16]基于該誤差函數(shù)分析了N和q的增加對相對誤差降低速度的影響,導(dǎo)出了在給定精度EPS時參數(shù)(N,q)優(yōu)化組合的自適應(yīng)選擇算法和計算流程,此處不再贅述。

3 時不變線性微分方程的精細積分方法

許多問題的數(shù)學(xué)模型都可采用時不變線性微分方程組來描述,即在一般性描述式(1)中,A是時不變矩陣,且f與狀態(tài)v不相關(guān),表示為

(25)

(k=0,1,2,…)

(26)

在區(qū)段[tk,tk+1]上利用Duhamel積分式(2),得到時不變線性微分方程式(25)的離散遞推列式為

(27)

式中

(28)

3.1 解析積分方法

鐘萬勰首先在文獻[2]中給出了當非齊次項f在區(qū)段[tk,tk+1]為線性函數(shù)或采用線性近似時的解析積分格式。隨后,林家浩等[17,18]通過引入特解vp(t)將式(27)表示為

vk+1=T[vk-vp(tk)]+vp(tk+1)

(29)

并進一步導(dǎo)出了在時間區(qū)段[tk,tk+1]內(nèi)非齊次項f(t)為多項式、指數(shù)函數(shù)、三角函數(shù)及其組合函數(shù)時對應(yīng)的特解vp(t)的解析格式,介紹如下。

(1) 當非齊次項f(t)為線性函數(shù),即

f(t)=f0+f1·(t-tk)

(30)

其對應(yīng)的特解為

vp(t)=(A-1+I·t)(-A-1f1)-A-1(f0-f1·tk)

(31)

將式(31)代入式(29),進行整理即得到如下解析積分格式

vk+1=T[vk+A-1(f0+A-1f1)]-

(32)

(2) 當非齊次項f(t)為三角函數(shù),即

f(t)=f1sin(ωt)+f2cos(ωt)

(33)

其對應(yīng)的特解為

vp(t)=asin(ωt)+bcos(ωt)

(34)

式中

(35)

(3) 當f(t)為多項式與三角函數(shù)的組合,即

f(t)=(f0+f1·t+f2·t2)×

[αsin(ωt)+βcos(ωt)]

(36)

其對應(yīng)的特解為

vp(t)=(a0+a1·t+a2·t2)sin(ωt)+(b0+b1·t+b2·t2)cos(ωt)

(37)

式中

(i=2,1,0)

(38)

(39)

(4) 當f(t)為指數(shù)與三角函數(shù)的組合,即

f(t)=eαt[f1sin(ωt)+f2cos(ωt)]

(40)

其對應(yīng)的特解為

vp(t)=eαt[asin(ωt)+bcos(ωt)]

(41)

式中

(42)

以上α,β,f0,f1,f2,a0,a1,a2,b0,b1和b2等在時間區(qū)段[tk,tk+1]內(nèi)均為常量。

當非齊次項具有上述函數(shù)形式時,解析積分格式本身沒有任何近似。對于一般形式的非齊次項,也可在時間區(qū)段[tk,tk+1]內(nèi)用上述函數(shù)進行近似,由于遞推列式簡單,得到了廣泛應(yīng)用。

值得指出的是,上述解析積分格式中含有矩陣的逆A-1或(ω2I+A2)-1,而求逆運算容易引起數(shù)值困難,特別是當矩陣奇異或接近奇異時。對此,文獻[19]指出對于奇異矩陣可先利用奇異值分解得到正交矩陣,再采用精細積分方法求解。文獻[20]進一步針對奇異Hamilton系統(tǒng)矩陣,利用Hamilton系統(tǒng)的共軛辛正交歸一關(guān)系將零本征值對應(yīng)子空間分離出來,再采用精細積分法求解。

此外,發(fā)展避免求逆的一般性算法也得到了學(xué)者們的關(guān)注,歸納起來主要有以下三類。

3.2 直接數(shù)值積分方法

為避免矩陣求逆運算導(dǎo)致的數(shù)值困難,許多學(xué)者將式(27)的Duhamel積分看作普通函數(shù)的積分,記

(43)

并分別采用常用的Simpson,Romberg,Cotes和Gauss等數(shù)值積分技術(shù)進行計算,然后與矩陣指數(shù)的精細積分相結(jié)合,構(gòu)造式(27)的離散遞推列式。

以采用Gauss數(shù)值積分為例[21,22],式(43)的數(shù)值積分如下

(44)

式中Ng為積分點的個數(shù),τi為Gauss積分點的坐標,wi是加權(quán)系數(shù)。將式(44)代入式(27)得到離散遞推列式為

(45)

式中

(46)

文獻[23]比較了多種直接數(shù)值積分的精度,指出Cotes和Gauss積分是保持精細積分算法高精度的較好方法,文獻[14,24]則進一步指出采用Padé逼近為基礎(chǔ)的矩陣指數(shù)精細方法是無條件穩(wěn)定的計算格式,可給出更高的精度。

通過直接數(shù)值積分與矩陣指數(shù)的精細積分方法相結(jié)合來構(gòu)造求解列式,不需要矩陣求逆運算,也無需對非齊次項進行數(shù)學(xué)擬合,具有很好的通用性。微分方程式(27)的求解精度主要取決于Duhamel數(shù)值積分的精度。同時,式(46)表明還需計算積分點上的矩陣指數(shù),增加了計算量。但這類方法沒有充分利用Duhamel積分函數(shù)中矩陣指數(shù)的特殊性。

3.3 增維精細積分方法

文獻[25]首先將增維技術(shù)應(yīng)用于非齊次微分方程的精細積分求解算法的構(gòu)造,指出如果非齊次項f為一組齊次線性常微分方程組的解,記為

(47)

(48)

式中

(49)

增維齊次方程(48)與原非齊次方程(25)是完全等價的,可直接采用矩陣指數(shù)的精細積分方法求解。

文獻[26]還提出了基于Taylor級數(shù)的增維精細積分方法(也稱為齊次擴容精細積分方法)。首先在時間區(qū)段[tk,tk+1]內(nèi)將非齊次項f進行Taylor級數(shù)展開,即設(shè)

f(t)=fm·tm+fm-1·tm-1+…+f1·t+f0

(50)

然后引入擴展狀態(tài)變量

(51)

(52)

此外,文獻[29]還給出了非齊次項為正/余弦函數(shù)形式的增維方法,記非齊次項為

f(t)=acos(ωt)+bsin(ωt)

(53)

式中a和b為不依賴時間的向量。引入擴展狀態(tài)變量

(54)

(55)

利用類似的思想,可以Legendre,Chebyshev,Hermite和Laguerre等正交多項式函數(shù)系[30,31]來構(gòu)造增維矩陣,并且文獻[31]還指出基于Legendre函數(shù)系構(gòu)造齊次方程具有更高的精度和收斂性。

增維精細積分方法將非齊次方程轉(zhuǎn)化為齊次方程,在實施精細積分的過程中不必進行矩陣求逆運算,同時保留了精細積分的高精度,其代價是增加了矩陣指數(shù)的計算量。

3.4 擴展精細積分方法

實際上,利用Duhamel積分中矩陣指數(shù)的特點,可以構(gòu)造出與解析積分格式一樣高效和精確的遞推算法,而且不需要矩陣求逆運算。即文獻[32]提出的擴展精細積分方法。

離散遞推列式(27)表明,其關(guān)鍵在于第二項Duhamel積分項的計算。首先,為了充分利用非齊次項本身的特點降低計算量,將非齊次項f(t)寫為

f(t)=Bfs(t)

(56)

式中 常系數(shù)矩陣B反映了非齊次項f(t)的特點,例如結(jié)構(gòu)動力學(xué)方程中作動力/干擾力的位置矩陣。由于fs(t)的維數(shù)小于f(t)的維數(shù),從而降低下面基函數(shù)響應(yīng)矩陣的維數(shù),提高計算效率。

將式(56)代入式(27),Duhamel積分項描述為

(57)

在區(qū)段[tk,tk+1]內(nèi)fs(t)采用特定基函數(shù)逼近為

(58)

(59)

(60)

顯然,基函數(shù)響應(yīng)矩陣僅與矩陣A的矩陣指數(shù)、基函數(shù)和區(qū)段長度有關(guān),體現(xiàn)了非齊次項在基函數(shù)模式下的響應(yīng)特征。對于時不變系統(tǒng),等長區(qū)段劃分的基函數(shù)的響應(yīng)矩陣只需計算一次。

將式(57,59)代入(27),即得到基于響應(yīng)矩陣的精細積分遞推公式為

(61)

利用基函數(shù)響應(yīng)矩陣表達式中含有矩陣指數(shù)的性質(zhì),文獻[32]也導(dǎo)出了基函數(shù)為多項式函數(shù)、指數(shù)函數(shù)、正/余弦函數(shù)及其組合函數(shù)形式時響應(yīng)矩陣的加法定理,同時給出了基函數(shù)響應(yīng)矩陣的精細積分方法,因此稱之為擴展精細積分方法。

與解析積分方法相比,擴展精細積分方法可同樣給出非齊次項為多項式函數(shù)、指數(shù)函數(shù)、正/余弦函數(shù)及組合函數(shù)形式時的精確解。更具優(yōu)勢的是,擴展精細積分方法不需要對矩陣求逆,具有更好的數(shù)值穩(wěn)定性和更廣泛的適用性,對進一步構(gòu)造非線性微分方程數(shù)值算法非常有意義[33,34]。

4 時變線性微分方程的精細積分方法

所謂時變線性微分方程組,就是常微分方程組的一般性描述式(1)的系統(tǒng)矩陣A隨時間變化,記為A(t),且f與狀態(tài)v不相關(guān),表示為

(62)

此時利用狀態(tài)傳遞矩陣和Duhamel積分表述的通解式(2)仍然成立,狀態(tài)傳遞矩陣也仍然滿足式(3～5)。在區(qū)段[tk,tk+1]上,遞推式表達為

(63)

本節(jié)主要講述以精細積分方法為基礎(chǔ)求解齊次時變系統(tǒng)對應(yīng)的狀態(tài)傳遞矩陣的計算方法,而非齊次線性時變微分方程的求解與非線性微分方程的求解類似,將在下節(jié)介紹。

式(62)對應(yīng)的齊次時變線性動力系統(tǒng)表示為

(64)

在區(qū)段[tk,tk+1]上,方程解的遞推式表達為

vk+1=Φ(tk+1,tk)vk

(65)

4.1 常值近似法

(66)

(67)

式中 定義李括號運算為

A1,A2?A1A2-A2A1

(68)

文獻[36]給出另一種4階Magnus級數(shù)近似結(jié)果

(69)

誤差理論分析表明,該方法的精度僅依賴于Magnus級數(shù)截斷的階次和區(qū)段長度。常值近似也是后面乘法攝動的基礎(chǔ)。

4.2 乘法攝動法

時變系統(tǒng)雖不能完全采用精細積分方法,但可采用攝動法進行求解,并且將定常系統(tǒng)的精細積分作為攝動的基礎(chǔ)。鐘萬勰等[37]首先提出了時變Hamilton系統(tǒng)的矩陣指數(shù)計算的保辛乘法攝動的思想,隨后譚述君等[38]進一步將其應(yīng)用于變系數(shù)Riccati方程的保辛攝動近似求解。在此基礎(chǔ)上,富明慧等[39]將乘法攝動思想做了進一步推廣,提出了一種高階乘法攝動方法,將時變系統(tǒng)狀態(tài)傳遞矩陣轉(zhuǎn)換為一系列矩陣指數(shù)的乘積,并可采用精細積分方法精確計算。

對齊次時變系統(tǒng)(64),取區(qū)段[tk,tk+1]上的一點tc,將時變矩陣A(t)分解為定常和時變兩部分,即

(70)

(71)

將式(71)代入式(64)可得到一階乘法攝動系統(tǒng),即

(72)

式中

(73)

這樣原系統(tǒng)(64)就轉(zhuǎn)化為一階攝動系統(tǒng)(72),其系數(shù)矩陣為關(guān)于(t-tc)的一階小量?？梢钥闯?當tc=tk時,式(73)即為文獻[38]提出的乘法攝動。

依次類推,繼續(xù)進行攝動。設(shè)攝動的最終次數(shù)為M,其狀態(tài)變換如下

(74)

相應(yīng)的第M階乘法攝動系統(tǒng)為

(75)

式中

(76)

將AM(t)分解為

(77)

(78)

式中α=tc-tk,β=tk+1-tc

(79)

文獻[39]還討論了tc=tk,tc=tk+1和tc=(tk+tk+1)/2時Φ(tk+1,tk)的精度和計算量,指出tc=(tk+tk+1)/2時算法性價比最高。

值得指出的是,如果系統(tǒng)矩陣A(t)是Hamilton矩陣,那么攝動后的系統(tǒng)仍然是Hamilton系統(tǒng),因此該方法成為一種高階保辛攝動方法,這對于Hamilton系統(tǒng)求解的保辛要求是非常重要的。

4.3 周期變系數(shù)系統(tǒng)

周期變系數(shù)系統(tǒng)是一般時變線性系統(tǒng)的特例,即式(64)的時變系統(tǒng)矩陣具有周期時變特性

A(t)=A(t+T)

(80)

式中T是周期的大小。

根據(jù)Floquet-Lyapunov理論,周期時變系統(tǒng)的穩(wěn)定性由Floquet轉(zhuǎn)移矩陣特征值的實部決定。Floquet轉(zhuǎn)移矩陣是狀態(tài)傳遞矩陣在一個周期端部的值,參照式(2)的描述,記為Φ(t0+T,t0)。Floquet轉(zhuǎn)移矩陣的計算,以往主要有兩類方法,一類是Hsu[40]提出的一系列將矩形波方法推廣于多維系統(tǒng)的近似方法;另一類是各種直接數(shù)值積分方法,其中,四階預(yù)估-校正Hamming法最為有效[41]。利用系數(shù)周期性變化的特點,文獻[42]將解析精細積分方法應(yīng)用于Floquet轉(zhuǎn)移矩陣的計算。然而該方法需要求逆運算,難以處理奇異或接近奇異矩陣情況,對此文獻[43]建立了相應(yīng)的擴展精細積分方法。

根據(jù)周期函數(shù)的Fourier級數(shù)展開式,式(80)描述的周期函數(shù)A(t)可采用s階Fourier級數(shù)近似

(81)

式中ωi=i2π/T,而A0與Di,Bi(i=1,2,…,s)均為常數(shù)矩陣。將式(81)代入式(64),得到

(82)

式中

(83)

將一個周期[t0,t0+T]的時間長度均勻劃分為M份,一系列等間距時刻為

(k=0,1,2,…,M)

(84)

在區(qū)段[tk,tk+1]內(nèi)對v(t)采用p階多項式近似

(85)

將式(85)代入式(83),得到

(86)

式中

(87)

在區(qū)段[tk,tk+1]上利用Duhamel積分,可以得到區(qū)段內(nèi)的遞推列式

(88)

(89)

式中

(90)

(91)

再利用傳遞矩陣性質(zhì)式(4),即可得到周期[t0,t0+T]內(nèi)任意時刻的狀態(tài)傳遞矩陣

(k=0,1,…,M-1)

(92)

式中Φ(t0,t0)=I。

文獻[43]分別給出了基于擴展精細積分的零階格式、一階格式和二階格式。這些格式與同階次的解析精細積分方法[42]相比,具有相同的精度和更高的效率,而且避免了矩陣求逆,具有更好的穩(wěn)定性。該方法也用于周期時變系統(tǒng)H2范數(shù)的計算[44]。

4.4 多項式變系數(shù)系統(tǒng)

多項式變系數(shù)系統(tǒng)是一般時變線性系統(tǒng)的另一種特例,在式(64)的時變系統(tǒng)矩陣可以采用多項式函數(shù)描述或近似表示為

(93)

式中Ai為時不變系數(shù)矩陣。文獻[45]引入新的變量和增維技術(shù),將其轉(zhuǎn)換為時不變系統(tǒng)的矩陣指數(shù),再采用精細積分方法精確求解。

首先,引入變換

τ=ρ(t-t0)/(tf-t0)

(94)

將原時變系統(tǒng)(64)轉(zhuǎn)換為如下單位時變系統(tǒng)

(95)

(i=1,2,…,m+r;m>r)

(96)

(97)

(98)

(99)

類似的,4.3節(jié)求解周期變系數(shù)系統(tǒng)狀態(tài)傳遞矩陣的擴展精細積分方法,也可以用于本節(jié)多項式變系數(shù)系統(tǒng)的狀態(tài)傳遞矩陣的計算。

5 非線性微分方程的精細積分方法

所謂非線性微分方程組,就是式(1)的右端含有狀態(tài)v的非線性項。不失一般性,可以將非線性項放在系統(tǒng)矩陣A中,記為A(v,t),f仍然可以表示成與狀態(tài)v不相關(guān)的形式表示為

(100)

式中 矩陣A(v,t)體現(xiàn)非線性和時變特征。顯然線性時變系統(tǒng)是其特例,因此,下面的數(shù)值算法也適用于含非齊次項的線性時變微分方程組的求解。

5.1 基于線性系統(tǒng)的精細積分方法

考慮到線性微分方程求解上的優(yōu)勢,一個自然的想法是將非線性部分的線性部分單獨提取出來,即將A(v,t)分解為

A(v,t)=A0+A1(v,t)

(101)

式中A0為定常系統(tǒng)矩陣,A1(v,t)表示非線性/時變部分,將其并入到非齊次項,將式(100)改寫為

(102)

式中

(103)

對于式(102)描述的非線性微分方程的解,在區(qū)段[tk,tk+1]內(nèi)采用Duhamel積分描述的遞推列式為

(104)

目前,對Duhamel積分中的非線性非齊次項在tk處通常采用多項式函數(shù)近似。文獻[46]給出了基于一次多項式近似并結(jié)合解析積分公式導(dǎo)出的精細積分算法。進一步,文獻[47]將其擴展到采用m次多項式近似的精細積分算法構(gòu)造。

(105)

(106)

多項式函數(shù)的Duhamel積分可以給出解析表達式,其遞推式如下

(107)

(108)

避免了解析積分格式(107)的矩陣求逆運算,但是以損失精度為代價。

實際上,第3節(jié)中面向時不變線性微分方程求解介紹的直接數(shù)值積分方法、增維精細積分方法和擴展精細積分方法等避免矩陣求逆運算的方法都可用于式(106)的求解。下面介紹文獻[32,33]基于擴展精細積分和外插法構(gòu)造的更為簡潔的非線性微分方程的精細算法。

基于擴展精細積分方法[32],當基函數(shù)取為多項式函數(shù)時,多項式響應(yīng)函數(shù)表示為

(j=0,1,…,m)

(109)

則式(106)的遞推公式變?yōu)?/p>

(110)

遞推式(110)唯一的近似是多項式近似,其與不同數(shù)值近似方法結(jié)合將導(dǎo)出不同的非線性微分方程數(shù)值算法。文獻[43,33]進一步導(dǎo)出了幾種單步法/多步法的顯式/隱式計算格式。以Adams線性多步法為例,如果利用tk-m,tk-m+1,…,tk-1,tk作為Lagrange插值節(jié)點,可導(dǎo)出基于擴展精細積分方法的顯式Adams的m步法計算格式,統(tǒng)一描述為

(111)

(112)

式中αj,l為顯式Adams的m步法中基函數(shù)響應(yīng)矩陣的組合系數(shù),表1給出了m=1,2,3步法時的組合系數(shù)。

利用tk-m+1,…,tk-1,tk,tk+1(插值點包含未知的tk+1)作為Lagrange插值節(jié)點,可導(dǎo)出基于擴展精細積分方法的隱式Adams的m步法,統(tǒng)一描述為

(113)

(114)

式中αj,l為隱式Adams的m步法中基函數(shù)響應(yīng)矩陣的組合系數(shù),表2給出了m=1,2,3步法時的組合系數(shù)。

表2 隱式Adams的m步法的αj,l

文獻[34,49]在擴展精細積分方法基礎(chǔ)上,將結(jié)構(gòu)動力學(xué)方程中的剛度陣項看作非齊次項,改寫為

(115)

(116)

(117)

顯然,利用上述結(jié)構(gòu)特點和分塊矩陣運算,矩陣指數(shù)和響應(yīng)矩陣只有一半的元素需要計算,同時狀態(tài)空間的遞推式(110)也可以利用分塊矩陣運算進一步簡化,從而降低計算量。

文獻[50]還采用四階Runge-Kutta法來計算式(104)右端的非齊次項的Duhamel積分,而采用精細積分方法來計算齊次部分,提出了精細Runge-Kutta法的計算列式為

(118)

式中

(119)

基于式(102)的算法,A0的選擇會影響到算法的精度和穩(wěn)定性。然而,如何選擇最優(yōu)的線性系統(tǒng)矩陣A0,尚未找到嚴格的方法。文獻[46]指出,為了后續(xù)算法的構(gòu)造,可以通過在A0中填充元素的方式使其滿秩。這似乎說明這些方法更適用于弱非線性問題,而對于強非線性問題,如何利用精細積分方法的優(yōu)點還是值得研究的問題。

5.2 基于同倫攝動的精細積分方法

同倫攝動方法是近年來一種求解非線性問題的漸近數(shù)值方法,具有固定的計算格式和推導(dǎo)方法。文獻[51]將精細積分技術(shù)和同倫攝動方法相結(jié)合,提出了基于精細積分技術(shù)的非線性動力學(xué)方程的同倫攝動法,具有較高的計算精度和效率。

針對非線性微分方程(102),構(gòu)造以下線性同倫函數(shù),即

v(t0)=v0

(120)

(121)

將式(121)代入式(120),并令方程ε同次冪的系數(shù)相等,得

(122)

顯然上述構(gòu)造的同倫函數(shù)的每一個攝動方程都是時不變線性微分方程,可用第3節(jié)介紹的精細積分方法求解。得到這一組攝動解之后,代入式(121)并令ε=1,即得到原非線性微分方程的解答。

5.3 增維為齊次時變系統(tǒng)的方法

(123)

則非線性微分方程(100)轉(zhuǎn)化為齊次方程形式,即

(124)

文獻[53]提出了一種在Minkowski空間中齊次化非線性系統(tǒng)的增維方法。對一般非線性系統(tǒng)

(125)

(126)

式中

(127)

(128)

式中G為一個Minkowski矩陣

(129)

這樣,原非線性動態(tài)系統(tǒng)就轉(zhuǎn)化為一個增維的Lie型動態(tài)系統(tǒng),可采用精細積分方法求解。文獻[54]則進一步采用Runge-Kutta Munthe-Kaas方法和精細積分方法構(gòu)造了簡潔的保群結(jié)構(gòu)積分算法。

6 大規(guī)模問題的精細積分方法

大規(guī)模工程問題的瞬態(tài)響應(yīng)分析,其相應(yīng)的系數(shù)矩陣通常具有很高的維度,若直接采用精細積分方法進行計算,需要計算高維度矩陣的相乘并破壞原始系數(shù)矩陣的稀疏性,從而導(dǎo)致消耗大量的計算時間和儲存空間。

6.1 子域精細積分法

為降低精細積分方法的計算量,鐘萬勰等[55-57]提出了子域精細積分方法。主要思想是利用系數(shù)矩陣具有窄帶寬的特點將結(jié)構(gòu)分為若干個子域,將計算轉(zhuǎn)化為一些子域上的精細積分計算。在此基礎(chǔ)上,陳飆松等[58]提出了子結(jié)構(gòu)精細積分方法,該方法將物理域劃分子結(jié)構(gòu),更適合于有限元法及并行計算?？紤]瞬態(tài)熱傳導(dǎo)問題的一階常微分方程組

(130)

式中C和K分別為熱容矩陣和熱傳導(dǎo)矩陣,T和F分別為節(jié)點溫度和熱載荷向量。假設(shè)結(jié)構(gòu)可以分為A和B兩個子結(jié)構(gòu),則式(130)可改寫為

(131)

令

(132)

則式(131)可以進一步改寫為

(133)

根據(jù)精細積分方法,式(133)的解可以表示為

(134)

在時間段[tk,tk+1]內(nèi),做如下近似

(p=A,B)

(135)

(136)

子域精細積分方法提出后,針對子域精細積分方法的改進和應(yīng)用做了大量工作。蔡志勤等[59]對子域精細積分方法的穩(wěn)定性進行了分析,證明了單點、兩點和三點子域精細積分及單點子域精細積分的蛙跳格式都是無條件穩(wěn)定的。曾文平[60,61]分別對二維和三維擴散方程提出了單點子域精細積分方法。金承日等[62,63]利用子域精細積分方法思想,針對對流方程初邊值問題,構(gòu)造出一組無條件穩(wěn)定的三層顯格式(蛙跳格式)和一組無條件穩(wěn)定的二層隱式格式,進而設(shè)計出求解該隱式格式的顯式迭代算法。賴永星等[64]基于單點子域精細積分思想,針對拋物線型熱傳導(dǎo)方程初邊值問題,提出多點子域積分的概念,推出一種多點子域積分的FTCS(Forward for Time,Center for Space)格式,并且說明了多點子域積分的FTCS格式具有比單點子域積分的FTCS格式收斂速度快的特點。Wu等[65]針對由相同結(jié)構(gòu)單元組成的周期結(jié)構(gòu)動力響應(yīng)問題,基于精細積分方法、子域格式和周期結(jié)構(gòu)的可重復(fù)性提出了一種子域精細積分方法,此方法特別適用于周期結(jié)構(gòu)動力響應(yīng)的求解。

6.2 Krylov子空間精細積分方法

Fung等[66]首先將Krylov子空間方法和精細積分方法相結(jié)合,提出了一種求解大規(guī)模瞬態(tài)問題的高效精細積分算法。對一個給定的n維方陣A以及非零向量v,定義m維Krylov子空間為

Km≡span{v,Av,A2v,…,Am-1v}

(137)

(1) 初始化

β=‖v‖2,v1=v/β

(138)

(2) 迭代過程

Doj=1,2,…,m-1

①w=Avj

② Doi=1,2,…,j

hi,j=wTvi

w=w-hi,jvi

③hj+1,i=‖w‖2,vj+1=w/hj+1,i

(139)

式中 向量e1的第一個元素為1,其余元素均為0。

在此基礎(chǔ)上,陳臻林[67]提出了一種結(jié)合Krylov子空間方法的精細時程積分法求解大型動力系統(tǒng),該方法利用增維技術(shù)避免了非齊次方程特解的求解,可處理任意形式載荷。Wang[68]在利用分段插值多項式模擬外載荷基礎(chǔ)上,提出了一種使用Ritz矢量的精細積分方法和一種改進的Krylov精細積分方法。Chen[69]針對具有任意激勵的二階微分方程,提出了一種改進的Krylov精細積分算法,該方法通過精確確定迭代子空間大小以解決復(fù)雜激勵問題。

6.3 基于物理特性的快速精細積分方法

高強等[70,71]根據(jù)結(jié)構(gòu)動力響應(yīng)的物理特點,從物理上說明了大規(guī)模動力系統(tǒng)對應(yīng)的小時間段矩陣指數(shù)應(yīng)該具有稀疏性,并基于該性質(zhì)提出了一種改進的快速精細積分方法。

空間離散后的結(jié)構(gòu)動力學(xué)方程為

(140)

式(140)可以寫為狀態(tài)空間形式,即

(141)

其中

(142)

(143)

式中

(144)

離散結(jié)構(gòu)動力學(xué)響應(yīng)的傳遞速度是有限的,下面利用這個性質(zhì)說明其對矩陣指數(shù)結(jié)構(gòu)的影響。根據(jù)式(143),如果外力為零,則式(143)可以改寫為

qk+1=T11qk+T12pk,pk+1=T21qk+T22pk

(145)

假設(shè)某個自由度上有擾動,在較小時間內(nèi)必然只能傳播到有限的自由度,而不是所有自由度。根據(jù)矩陣T11,T12,T21和T22的物理含義,則其一定是稀疏矩陣?；谝陨戏治?該方法對精細積分方法計算過程進行稀疏化過濾處理,即在循環(huán)計算矩陣指數(shù)T的增量部分時,在每一循環(huán)結(jié)束前都對矩陣指數(shù)內(nèi)由于計算誤差產(chǎn)生的小于容差的數(shù)進行過濾(即其值設(shè)為0),來提高矩陣指數(shù)的稀疏性,從而提高計算效率。根據(jù)類似思想,吳鋒等[72]針對雙曲熱傳導(dǎo)問題的物理特點,提出了一種求解大規(guī)模雙曲熱傳導(dǎo)問題的快速精細積分方法。

6.4 基于切比雪夫展開的高效率精細積分方法

考慮如下一類特殊形式的一階常微分方程組

(146)

式中 矩陣C為對稱正定稀疏矩陣,矩陣K為對稱半正定稀疏矩陣。如有限元離散化后的瞬態(tài)熱傳導(dǎo)問題的一階常微分方程組即滿足式(146)。

針對于這類特殊的常微分方程組問題,高強等[73]提出了切比雪夫展開高效數(shù)值方法。該方法基于矩陣指數(shù)的切比雪夫矩陣多項式展開,將矩陣指數(shù)-向量的乘積轉(zhuǎn)化為計算一系列切比雪夫矩陣多項式與向量的乘積,從而顯著提高了計算效率。

對矩陣C進行Cholesky分解C=LLT,則式(146)可改寫為

(147)

(148)

然后采用以下切比雪夫展開法計算式(148)的矩陣指數(shù)-向量乘積

(149)

式中

(150)

Anorm=(2A-λmaxI)/λmax

(151)

式中δn,0滿足當n=0時δn,0=1,而當n>0時δn,0=0;λmax為矩陣A的最大特征值;In(x)為n階第一類修正貝塞爾函數(shù);Tn(X)為n階第一類切比雪夫矩陣多項式,可按以下遞歸過程計算

(152)

(153)

式中

(154)

由式(154)可知,由于向量un和n有關(guān),計算式(153)右端級數(shù)的每一項都需要計算一次向量積分,根據(jù)實際問題中外載荷隨時間變化特點,文獻[73]給出了高效計算式(153)的方法。

穩(wěn)定性分析表明,對于任意時間步長,該方法近似算子的譜半徑均小于1,因此該方法具有無條件穩(wěn)定的特性。同時,其計算過程中只涉及矩陣-向量乘積計算并只需儲存少數(shù)向量,因此針對大規(guī)模瞬態(tài)熱傳導(dǎo)問題,該方法在計算效率和內(nèi)存消耗上都具有很大的優(yōu)勢。

6.5 周期結(jié)構(gòu)的高效率精細積分方法

基于周期結(jié)構(gòu)動力響應(yīng)的物理特點,高強等[74-76]提出了周期結(jié)構(gòu)動力響應(yīng)分析的高效率精細積分方法。

假設(shè)周期結(jié)構(gòu)的一個單胞自由度數(shù)為d,并將矩陣指數(shù)T按照方程(144)分塊,則周期結(jié)構(gòu)對應(yīng)的矩陣指數(shù)具有以下特點

(155)

由式(155)可知,周期結(jié)構(gòu)對應(yīng)的矩陣指數(shù)中有很多相同的矩陣元素,為提高計算效率,這些相同的元素不應(yīng)該全部計算和存儲。另一個影響因素是能量在結(jié)構(gòu)中的傳遞速度是有限的。假設(shè)某個自由度上有擾動,在較小步長內(nèi),必然只能傳播到有限的自由度,而不可能傳播到所有自由度。

綜合以上兩個性質(zhì),則矩陣塊T11,T12,T21和T22都具有如下的帶狀結(jié)構(gòu),即

[]

(156)

式中 標量ai(i=1,2,…,d)對應(yīng)于矩陣的對角線元素,向量bi(i=1,2,…,d)的第一個元素不為零,而向量ci(i=1,2,…,d)的最后一個元素不為零。因此,要計算周期結(jié)構(gòu)矩陣指數(shù),只需要計算和存儲ai,bi和ci(i=1,2,…,d)即可。

下面以圖1(a)所示的含有N個單胞的一維周期結(jié)構(gòu)來說明ai,bi和ci(i=1,2,…,d)的計算方法。假設(shè)在一個積分步長內(nèi),能量傳播的速度不超過m個單胞,則可構(gòu)造如圖1(b)所示的由2m+1個單胞組成的代表單胞結(jié)構(gòu)。這兩個結(jié)構(gòu)的邊界都不會對第m+1個單胞上的位移和動量響應(yīng)產(chǎn)生影響,所以兩個結(jié)構(gòu)對應(yīng)節(jié)點的響應(yīng)相同,而結(jié)構(gòu)(a)中第2m+2到第N個單胞對應(yīng)的響應(yīng)為零。

圖1

因此,基于以上周期結(jié)構(gòu)動力響應(yīng)矩陣指數(shù)的代數(shù)結(jié)構(gòu)分析,計算結(jié)構(gòu)(b)的矩陣指數(shù),并從恰當?shù)奈恢锰崛?shù)據(jù)即可得到式(156)的ai,bi和ci(i=1,2,…,d)。從而大大降低了矩陣指數(shù)的計算量和儲存量。

基于類似思想,高強等[77,78]提出了求解一維和二維周期結(jié)構(gòu)瞬態(tài)熱傳導(dǎo)問題的高效率數(shù)值方法。在此基礎(chǔ)上,高強等[79,80]針對含有移動熱源的周期結(jié)構(gòu)和缺陷周期結(jié)構(gòu)瞬態(tài)熱傳導(dǎo)問題提出了高效數(shù)值方法。

7 精細積分思想的擴展應(yīng)用

精細積分方法的要點在于2N類算法和增量存儲?；谶@兩個要點,精細積分方法也擴展用于許多其他問題,解決了傳統(tǒng)方法在計算這些問題時精度低的困境。

7.1 兩點邊值問題

基于計算結(jié)構(gòu)力學(xué)和最優(yōu)控制之間的模擬理論,鐘萬勰等[3,81,82]進一步將精細積分方法擴展至兩點邊值問題。

設(shè)整個兩端邊值問題的求解域為[za,zb],兩端邊值問題的一階2n維齊次方程可描述為

dv/dz=Av,v=[qT,pT]T

(157)

式中A為2n維的系統(tǒng)矩陣,v為2n維的待求狀態(tài)向量。在邊界za和zb處各有n個已知量,不失一般性,本文設(shè)為qa=q(za)和pb=p(zb)。對于兩端邊值問題精細積分法雖然不能直接使用,但是其增量存儲和基于加法定理的2N類算法的思想仍然適用。式(157)通?？梢圆捎镁仃囍笖?shù)方法或區(qū)段混合能方法進行求解。

(1) 矩陣指數(shù)方法

設(shè)兩端邊值問題的求解域為[za,zb],即求解域長度為L=zb-za,則方程(157)的解可表示為

(158)

式中Φ(L)=eAL可采用精細積分求解。將Φ(L)寫成分塊矩陣形式,則有

(159)

式(159)將兩端的狀態(tài)向量聯(lián)系起來。由于已知qa和pb,則重新整理式(159)可得

(160)

(2) 區(qū)段混合能方法

區(qū)段混合能方法常用于最優(yōu)控制[84],鐘萬勰[83]提出了用區(qū)段混合能求解非對稱Riccati方程的方法,這一方法也可用到求解兩端邊值問題(157)。

考慮到力學(xué)問題中,A為Hamilton矩陣的情況比較常見,本文以Hamilton矩陣為例進行介紹。將式(157)寫成分塊矩陣形式,即

(161)

式中B和D為對稱矩陣。此時式(160)具有如下形式

(162)

式中F,G和Q是區(qū)段長度L=zb-za的矩陣函數(shù)。參考精細積分方法的思想,首先將求解區(qū)域等分為2N個小區(qū)段,在每個小區(qū)段τ=2-NL上,F,G和Q可用M項Taylor級數(shù)近似展開為

(163)

式中θk,γk和φk為Taylor展開的系數(shù)矩陣

θ1=B,γ1=D,φ1=C

(164)

而k=2,3,…,M時

(165)

同時,還存在增量加法定理

(166)

兩點邊值問題的精細積分方法提出后,首先在LQ控制[89,90]、Riccati方程[82]、卡爾曼-布西濾波[91]和H∞控制系統(tǒng)[92,93]等方面得到廣泛應(yīng)用,文獻[84,94]等對這些應(yīng)用作了詳細介紹。兩點邊值問題的精細積分方法還廣泛用于波傳播問題[95-100]、分層地基動力響應(yīng)問題[101,102]、電磁波反射和投射問題[103]和奇異攝動邊值問題[104,105]等方面。

7.2 橢圓函數(shù)

橢圓函數(shù)是由橢圓積分定義的一類函數(shù),廣泛應(yīng)用于工程問題。橢圓函數(shù)的直接計算一般比較困難,但可進行冪級數(shù)展開,且也有加法定理的重要性質(zhì),因此可采用精細積分方法思想進行求解[106]。

以Jacobi橢圓函數(shù)為例,介紹橢圓函數(shù)的精細積分方法。Jacobi橢圓函數(shù)的定義可從第二類Legendre橢圓積分引出,該積分為

(167)

式中 參數(shù)m0稱為模數(shù)。Jacobi橢圓函數(shù)sn(u,m0)=z定義為上述第二類Legendre橢圓積分的反函數(shù)。而另外兩個基本Jacobi橢圓函數(shù)可分別定義為

(168)

式中 省略了參數(shù)m0,這是文獻常用寫法。

這三個Jacobi橢圓函數(shù)有如下的冪級數(shù)展開

(169)

和如下加法定理

(170)

如果計算Jacobi橢圓函數(shù)在u點的值,可令τ=u/2N,此時τ將為非常小的數(shù)值,采用較少項數(shù)的冪級數(shù)展開即可較為精確地近似

(171)

再根據(jù)式(170,171)可得

(172)

利用式(172)迭代N次,可以得到

sn(u)=SN, cn(u)=1+CN, dn(u)=1+DN

(173)

實際的數(shù)值實驗表明,精細積分方法對于m0取任意值均可給出極高精度的計算結(jié)果。目前,橢圓函數(shù)的精細積分方法已成功應(yīng)用在海床動態(tài)響應(yīng)[107]、Schwarz-Christoffel變換模型求解[108]和淺水波方程求解[109]等問題上。

7.3 矩陣正/余弦

矩陣正/余弦也是工程中經(jīng)常用到的一種矩陣函數(shù)。對于較小規(guī)模問題,一般可以采用Jordan標準型及級數(shù)展開法精確地求解。但當矩陣維度較大時,Jordan方法并不合適。精細積分法也可用于矩陣正/余弦的計算,并具有十分出色的表現(xiàn)[110]。

取適當?shù)恼麛?shù)N,則矩陣A的正/余弦化為

sinA=sin(2NB), cosA=cos(2NB)

(174)

式中B=A/2N。記

sin(2iB)=Si, cos(2iB)=I+Ci

(175)

因B的數(shù)值非常小,因此根據(jù)Taylor級數(shù)展開有

sinB≈B-B3/3!+B5/5!=S0

cosB≈I-B2/2!+B4/4!=I+C0

(176)

式中I為單位矩陣。然后利用矩陣正/余弦的2倍角公式,可得到如下的遞推關(guān)系

(177)

通過式(177)對增量進行N次的計算后,最終可獲得

sinA=SN, cosA=I+CN

(178)

此時計算得到的增量CN已經(jīng)不再是一個很小的量,當其與單位矩陣相加時,計算機的舍入誤差已經(jīng)可以忽略不計,能夠保證計算的精度。

矩陣正/余弦的精細積分方法無需計算矩陣的特征值和特征向量,而且對實復(fù)數(shù)矩陣均適用,并可以給出計算機意義上的精確解,非常適合大規(guī)模矩陣的計算。

7.4 分數(shù)階微分方程

一些復(fù)雜力學(xué)物理過程往往具有明顯的記憶、遺傳和路徑依賴性質(zhì)[111],由于分數(shù)階微積分具有良好的時間記憶性[112],已經(jīng)廣泛應(yīng)用于這些復(fù)雜力學(xué)物理過程中[113,114]。然而,求解分數(shù)階微分方程十分困難,常用的方法包括多項式法[115,116]、模型降階法[117]和龍格-庫塔法[118]等,但這些方法的精度均不夠理想。眾多學(xué)者也嘗試將精細積分方法用于分數(shù)階微分方程,以期提高求解精度。一類是將分數(shù)階微分動力方程近似轉(zhuǎn)化為整數(shù)階微分方程,之后再利用精細積分方法求解對應(yīng)的近似整數(shù)階微分方程[119,120];另一類是直接對分數(shù)階微分方程構(gòu)建精細積分方法[121]。本節(jié)對第二類研究進行介紹。

分數(shù)階常微分方程的一般形式為

D(α)x(t)=Ax(t)

(0<α<1)

(179)

式中x(t)=[x1(t),x2(t),…,xn(t)]T,D為導(dǎo)數(shù)算子,A∈Rn×n。

利用拉氏變換及逆變換可得到式(179)的解為

x(t)=Eα,1(Atα)x(0)

(180)

式中Eα,β(X)為雙變量Mittag-Leffler函數(shù)(簡稱ML函數(shù))

(181)

式中α和β為兩個變量,Γ(·)表示伽馬函數(shù)。記

(182)

式中

ak=1/Γ(1+kα)

(k=0,1,2,…)

(183)

(184)

于是有

(185)

式中

(186)

顯然ML函數(shù)不滿足加法定理,因此建立F(2X)和F(X)的關(guān)系時,需要引入修正項R(X)。記

(187)

將式(184,186)代入式(187)即可給出修正項函數(shù)R(X),顯然其可以表示為

(188)

式中 系數(shù)bk可用ak(k=2,3,…,∞)表示。實際數(shù)值計算時,函數(shù)R(X)的表達式可近似地取為有限項。當Q1和Q2的多項式定義式分別截取m1和m2項,所得函數(shù)R(X)的多項式次數(shù)為max(m2,2m1),為提高精度,一般取m2≥m1。

于是可得到如下遞推關(guān)系式(i=1,2,…)

F(2iX)=I+Qi+1(2iX)

(189)

(190)

F(Aτ)=I+Q1(Aτ)

(191)

Q1(Aτ)≈Aτ/Γ(1+α)+(Aτ)2/Γ(1+2α)+(Aτ)3/Γ(1+3α)+…+(Aτ)s/Γ(1+sα)

(192)

式中s為整數(shù),類比精細積分方法,可取s=[4/α],[·]表示Gauss函數(shù),Q1是一個小量矩陣。

然后應(yīng)用遞推關(guān)系式(189,190),最終可得到ML方程(181)的精細積分方法的求解表達式,即

(193)

再結(jié)合式(180,193)即可解分數(shù)階常微分方程(179)。

由于分數(shù)階常微分方程解析解中的ML函數(shù)與常規(guī)的矩陣指數(shù)不同,不滿足加法定理,因此分數(shù)階微分方程的精細積分方法需要在迭代關(guān)系中加入修正項,同時因為修正項的存在,導(dǎo)致分數(shù)階微分方程的計算量增加,其計算量為普通整數(shù)階精細積分方法的p倍,其中p為多項式R(X)的次數(shù)。

7.5 病態(tài)代數(shù)方程

富明慧等[122]將病態(tài)代數(shù)方程歸結(jié)于一個常微分方程初值問題的極限形式,并結(jié)合精細積分方法,提出了一種基于精細積分思想的求解病態(tài)方程組的高效方法。

對于正定的系數(shù)矩陣A,由于

(194)

定義函數(shù)

(195)

因為

[I+e-Aτ]F(τ)

(196)

重復(fù)上述過程k次,可得

(197)

由方程(197)可知,隨著迭代次數(shù)的不斷增加,積分的上限以指數(shù)形式增加,因此該迭代方程可以高效地逼近逆陣A-1。實際計算時,迭代的初始值τ取為一非常小的數(shù)值,則F(τ)可用Taylor展開的前幾項近似計算,如可取

F(τ)≈Iτ-Aτ2/2+A2τ3/6-A3τ4/24

(198)

同時e-Aτ也只保留Taylor展開前幾項進行近似計算

e-Aτ≈I-Aτ+(Aτ)2/2-(Aτ)3/6=I+T0

(199)

將式(198,199)代入式(197),并注意對s=1,2,…有如下迭代格式,

(200)

即可得到病態(tài)系數(shù)矩陣A的逆A-1的高精度解。

理論上,隨著迭代次數(shù)的不斷增加,精細積分法的結(jié)果應(yīng)該不斷逼近理論解。然而矩陣乘積引起的計算誤差會隨著迭代的進行而積累,進而導(dǎo)致精度迅速下降甚至出現(xiàn)溢出。為此,富明慧等[123]針對該問題提出了相應(yīng)的迭代收斂準則。郝強等[124]引入了主元加權(quán)的思想,提高了精細積分的計算精度。王慧蓉等[125]將病態(tài)系數(shù)矩陣A進行分裂,以減小病態(tài)系數(shù)矩陣的條件數(shù)。此后,一些學(xué)者又對病態(tài)代數(shù)方程的精細積分方法的計算效率進行了一定程度的優(yōu)化,如張文志等[126]改進了迭代的格式,進一步降低了采用精細積分方法求解病態(tài)代數(shù)方程的計算量;富明慧等[127]利用范數(shù)均衡預(yù)處理法對病態(tài)系數(shù)矩陣A進行預(yù)處理,提升了精細積分方法計算效率。

8 結(jié) 論

矩陣指數(shù)的精細積分方法可以給出計算機意義上的精確解,得益于2N類算法和增量存儲兩個核心要點。2N類算法使得精細劃分和高效計算成為可能,而增量存儲則保證在執(zhí)行2N類算法的合并過程中避免計算機舍入誤差成為影響精度的主要因素,這也是矩陣指數(shù)精細積分方法成功的關(guān)鍵。

矩陣指數(shù)在微分方程的求解中具有重要的地位,因此在矩陣指數(shù)的精細積分方法取得成功之后,迅速用于線性/非線性微分方程數(shù)值算法的構(gòu)造,得到了一系列基于精細積分的高精度、高穩(wěn)定性數(shù)值算法,極大地豐富了微分方程的數(shù)值計算方法。同時,精細積分方法的思想也擴展應(yīng)用于很多其他問題,包括兩點邊值問題、橢圓函數(shù)、矩陣正/余弦函數(shù)、病態(tài)代數(shù)方程以及分數(shù)階微分方程等問題高精度求解算法的構(gòu)造中。這些已有工作很好展現(xiàn)出精細積分方法的特色和優(yōu)勢。

進一步,精細積分方法可望在下面幾個方向得到發(fā)展和突破。 (1) 保辛數(shù)值離散已成為Hamilton動力系統(tǒng)數(shù)值算法設(shè)計的重要原則,結(jié)合精細積分方法,可望設(shè)計出性能更高的保辛算法。 (2) 精細積分方法在弱非線性問題求解中表現(xiàn)出了優(yōu)異的性能,如何將精細積分方法應(yīng)用于強非線性問題,并構(gòu)造高性能的求解算法是一個挑戰(zhàn)。 (3) 復(fù)雜大規(guī)模問題的高效分析是解決工程問題的關(guān)鍵,在保證精細積分方法精度優(yōu)勢的同時提高計算效率是一個值得關(guān)注的問題。

致謝：謹以此文紀念鐘萬勰先生九十誕辰!