姚偉岸, 高 強(qiáng), 譚述君, 吳 鋒

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

1 引 言

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

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

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

2 精細(xì)積分方法的基本思想

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

(1)

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

(2)

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

(1)Φ(t,t)=I

(3)

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

(4)

(5)

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

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

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

(6)

式中 矩陣指數(shù)

(7)

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

2.1 精細(xì)積分方法的提出

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

(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陣是一個(gè)小量的矩陣,為

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

(11)

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

(12)

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

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

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

(13)

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

(14)

T=I+Ta

(15)

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

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

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

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

(16)

式中

(17)

(18)

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

參照精細(xì)積分方法,將精細(xì)區(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í)行過程中只存儲(chǔ)和計(jì)算增量部分,就得到了基于Padé級(jí)數(shù)近似的矩陣指數(shù)精細(xì)積分方法。

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

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

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

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

(23)

式中 EPS為指定精度

(24)

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

3 時(shí)不變線性微分方程的精細(xì)積分方法

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

(25)

(k=0,1,2,…)

(26)

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

(27)

式中

(28)

3.1 解析積分方法

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

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

(29)

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

(1) 當(dāng)非齊次項(xiàng)f(t)為線性函數(shù),即

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

(30)

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

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

(31)

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

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

(32)

(2) 當(dāng)非齊次項(xiàng)f(t)為三角函數(shù),即

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

(33)

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

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

(34)

式中

(35)

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

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

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

(36)

其對(duì)應(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) 當(dāng)f(t)為指數(shù)與三角函數(shù)的組合,即

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

(40)

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

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

(41)

式中

(42)

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

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

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

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

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

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

(43)

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

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

(44)

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

(45)

式中

(46)

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

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

3.3 增維精細(xì)積分方法

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

(47)

(48)

式中

(49)

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

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

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

(50)

然后引入擴(kuò)展?fàn)顟B(tài)變量

(51)

(52)

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

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

(53)

式中a和b為不依賴時(shí)間的向量。引入擴(kuò)展?fàn)顟B(tài)變量

(54)

(55)

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

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

3.4 擴(kuò)展精細(xì)積分方法

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

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

f(t)=Bfs(t)

(56)

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

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

(57)

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

(58)

(59)

(60)

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

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

(61)

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

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

4 時(shí)變線性微分方程的精細(xì)積分方法

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

(62)

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

(63)

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

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

(64)

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

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

(65)

4.1 常值近似法

(66)

(67)

式中 定義李括號(hào)運(yùn)算為

A1,A2?A1A2-A2A1

(68)

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

(69)

誤差理論分析表明,該方法的精度僅依賴于Magnus級(jí)數(shù)截?cái)嗟碾A次和區(qū)段長(zhǎng)度。常值近似也是后面乘法攝動(dòng)的基礎(chǔ)。

4.2 乘法攝動(dòng)法

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

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

(70)

(71)

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

(72)

式中

(73)

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

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

(74)

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

(75)

式中

(76)

將AM(t)分解為

(77)

(78)

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

(79)

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

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

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

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

A(t)=A(t+T)

(80)

式中T是周期的大小。

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

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

(81)

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

(82)

式中

(83)

將一個(gè)周期[t0,t0+T]的時(shí)間長(zhǎng)度均勻劃分為M份,一系列等間距時(shí)刻為

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

(84)

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

(85)

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

(86)

式中

(87)

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

(88)

(89)

式中

(90)

(91)

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

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

(92)

式中Φ(t0,t0)=I。

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

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

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

(93)

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

首先,引入變換

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

(94)

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

(95)

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

(96)

(97)

(98)

(99)

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

5 非線性微分方程的精細(xì)積分方法

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

(100)

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

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

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

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

(101)

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

(102)

式中

(103)

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

(104)

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

(105)

(106)

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

(107)

(108)

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

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

基于擴(kuò)展精細(xì)積分方法[32],當(dāng)基函數(shù)取為多項(xiàng)式函數(shù)時(shí),多項(xiàng)式響應(yīng)函數(shù)表示為

(j=0,1,…,m)

(109)

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

(110)

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

(111)

(112)

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

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

(113)

(114)

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

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

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

(115)

(116)

(117)

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

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

(118)

式中

(119)

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

5.2 基于同倫攝動(dòng)的精細(xì)積分方法

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

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

v(t0)=v0

(120)

(121)

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

(122)

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

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

(123)

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

(124)

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

(125)

(126)

式中

(127)

(128)

式中G為一個(gè)Minkowski矩陣

(129)

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

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

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

6.1 子域精細(xì)積分法

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

(130)

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

(131)

令

(132)

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

(133)

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

(134)

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

(p=A,B)

(135)

(136)

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

6.2 Krylov子空間精細(xì)積分方法

Fung等[66]首先將Krylov子空間方法和精細(xì)積分方法相結(jié)合,提出了一種求解大規(guī)模瞬態(tài)問題的高效精細(xì)積分算法。對(duì)一個(gè)給定的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的第一個(gè)元素為1,其余元素均為0。

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

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

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

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

(140)

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

(141)

其中

(142)

(143)

式中

(144)

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

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

(145)

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

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

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

(146)

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

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

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

(147)

(148)

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

(149)

式中

(150)

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

(151)

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

(152)

(153)

式中

(154)

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

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

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

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

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

(155)

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

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

[]

(156)

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

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

圖1

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

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

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

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

7.1 兩點(diǎn)邊值問題

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

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

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

(157)

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

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

設(shè)兩端邊值問題的求解域?yàn)閇za,zb],即求解域長(zhǎng)度為L(zhǎng)=zb-za,則方程(157)的解可表示為

(158)

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

(159)

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

(160)

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

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

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

(161)

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

(162)

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

(163)

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

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

(164)

而k=2,3,…,M時(shí)

(165)

同時(shí),還存在增量加法定理

(166)

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

7.2 橢圓函數(shù)

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

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

(167)

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

(168)

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

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

(169)

和如下加法定理

(170)

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

(171)

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

(172)

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

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

(173)

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

7.3 矩陣正/余弦

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

取適當(dāng)?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級(jí)數(shù)展開有

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

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

(176)

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

(177)

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

sinA=SN, cosA=I+CN

(178)

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

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

7.4 分?jǐn)?shù)階微分方程

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

分?jǐn)?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ù)(簡(jiǎn)稱ML函數(shù))

(181)

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

(182)

式中

ak=1/Γ(1+kα)

(k=0,1,2,…)

(183)

(184)

于是有

(185)

式中

(186)

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

(187)

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

(188)

式中 系數(shù)bk可用ak(k=2,3,…,∞)表示。實(shí)際數(shù)值計(jì)算時(shí),函數(shù)R(X)的表達(dá)式可近似地取為有限項(xiàng)。當(dāng)Q1和Q2的多項(xiàng)式定義式分別截取m1和m2項(xiàng),所得函數(shù)R(X)的多項(xiàng)式次數(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ù),類比精細(xì)積分方法,可取s=[4/α],[·]表示Gauss函數(shù),Q1是一個(gè)小量矩陣。

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

(193)

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

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

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

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

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

(194)

定義函數(shù)

(195)

因?yàn)?/p>

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

(196)

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

(197)

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

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

(198)

同時(shí)e-Aτ也只保留Taylor展開前幾項(xiàng)進(jìn)行近似計(jì)算

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

(199)

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

(200)

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

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

8 結(jié) 論

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

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

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

致謝：謹(jǐn)以此文紀(jì)念鐘萬勰先生九十誕辰!