1、大規(guī)模動(dòng)力系統(tǒng)改進(jìn)的快速精細(xì)積分方法高強(qiáng)，吳鋒，張洪武林家浩，鐘萬(wàn)勰(大連理工大學(xué)，工業(yè)裝備結(jié)構(gòu)分析國(guó)家重點(diǎn)實(shí)驗(yàn)室，工程力學(xué)系，遼寧大連116024)摘要：木文提出一種針對(duì)大規(guī)模動(dòng)力系統(tǒng)的改進(jìn)的快速精細(xì)積分方法(fpim)。以精細(xì)積分方 法為基礎(chǔ)，利用大規(guī)模動(dòng)力系統(tǒng)矩陣的稀疏性和動(dòng)力問(wèn)題的物理特性，分析了矩陣指數(shù)的特 殊結(jié)構(gòu)，并基于此給出一種計(jì)算大規(guī)模動(dòng)力系統(tǒng)矩陣指數(shù)及其動(dòng)力響應(yīng)的高效率方法。關(guān)鍵詞：動(dòng)力系統(tǒng)；稀疏矩陣；精細(xì)積分；矩陣指數(shù)；快速算法1引言對(duì)于線(xiàn)彈性結(jié)構(gòu)的動(dòng)力學(xué)問(wèn)題，比較成熟和常用的時(shí)域積分方法是newmark 方法111和runge-kutta方法|2_4,由于計(jì)算穩(wěn)定性和精度

2、的要求，這兩種方法的積 分步長(zhǎng)一般都取得比較小，計(jì)算量比較大。鐘萬(wàn)勰提出了精細(xì)積分方法ls7j，這 種方法計(jì)算精度非常高，穩(wěn)定性好，允許采用很大的積分步長(zhǎng)，特別是在處理剛 性問(wèn)題時(shí)具有明顯優(yōu)勢(shì)。精細(xì)積分方法提出后，得到了很多發(fā)展但是這種 方法的一個(gè)缺點(diǎn)是在處理規(guī)模很大的系統(tǒng)時(shí)，由于計(jì)算矩陣指數(shù)的計(jì)算量比較 大，效率是一個(gè)主要問(wèn)題。本文針對(duì)大規(guī)模動(dòng)力系統(tǒng)提出一種計(jì)算其動(dòng)力響應(yīng)的高效率算法。本文以精 細(xì)積分方法為基礎(chǔ)15_71,利用動(dòng)力問(wèn)題的物理特性，利用一個(gè)關(guān)鍵思想，也就是 結(jié)構(gòu)動(dòng)力問(wèn)題能量傳遞的有限性，來(lái)降低矩陣指數(shù)的計(jì)算量。利用上述思想，本 文分析了大規(guī)模動(dòng)力結(jié)構(gòu)對(duì)應(yīng)矩陣指數(shù)的稀疏性質(zhì)，并棊

3、于此給出一種計(jì)算矩陣 指數(shù)的高效率方法。在高效和精確計(jì)算矩陣指數(shù)的基礎(chǔ)上，給出了人規(guī)模動(dòng)力系 統(tǒng)響應(yīng)的高效率和高精度算法。2動(dòng)力系統(tǒng)的精細(xì)積分方法假設(shè)系統(tǒng)的剛度矩陣、阻尼矩陣和質(zhì)量矩陣分別為k, c和m,那么結(jié)構(gòu) 動(dòng)力學(xué)方程為a 然科學(xué)基金(10902020, 10632030, 10721062, 2009aa04450i),遼寧省博士啟動(dòng)基金(20081091),遼寧省重 點(diǎn)實(shí)驗(yàn)室項(xiàng)h(2009s018)，鐵道部科技研宂開(kāi)發(fā)課題(2010t001-c)，大連理工大學(xué)青年教師培養(yǎng)拈金。通訊作者：張洪成，人連理工人學(xué)工程力學(xué)系，電i5: 86 411 84706249; e-mail addr

4、ess: zhanghw12a22a '(ie o其中f(r)為外力。方程(1)可以寫(xiě)為狀態(tài)空間形式，即(2)hv +其中(3)其中p=為動(dòng)量。數(shù)值積分時(shí)，將積分區(qū)間等分成步長(zhǎng)為a的積分區(qū)間，即r0 = 0，t= h，l，tk= kh，tk+=(+ l)/z，l(4)若記v=v(，j(5)則方程(2)的解可以表示為v=tv d:exp鍵g- x)f(tk + x)d.x(6)其中t = exp(h/z)(7)通過(guò)方程(6)計(jì)算結(jié)構(gòu)的動(dòng)力響應(yīng)，需要解決兩個(gè)主要問(wèn)題，一是要準(zhǔn)確高 效的計(jì)算方程(7)定義的矩陣指數(shù)，二是要計(jì)算方程(6)中的積分。對(duì)于常見(jiàn)的外 載荷，方程(6)中積分大多可解析

5、積分，例如(1)如果外力在一個(gè)積分時(shí)間步長(zhǎng)內(nèi)為多項(xiàng)式，即f(z)= f()+f,z+l+f,/w，則tv, + tofo 4-t,f, 4-l 4-mm其屮，(=0, 1，l，m)cti2- t22),t(;2 = mti2(9)平卜膽匕，(k= 1, 2, l , m)上式中丁|7和1；7分別為矩陣指數(shù)t對(duì)應(yīng)的分塊矩陣，即t12(2)如果外力在一個(gè)積分時(shí)間步長(zhǎng)內(nèi)為簡(jiǎn)諧變化，即f (z)= rexp(/vvr)，則v,+ l=tv,+ exp(zw,x-i -捫 1 篆二)! - i: r 因此，上述計(jì)算過(guò)程的關(guān)鍵是計(jì)算矩陣指數(shù)。目前，計(jì)算矩陣指數(shù)最好的方 法是精細(xì)積分方法|5"71

6、。但是，精細(xì)積分法計(jì)算矩陣指數(shù)的計(jì)算量比較大，即使 對(duì)于h為稀疏矩陣的情況，也很難利用其稀疏性，得到的矩陣指數(shù)可能是滿(mǎn)陣。3改進(jìn)的快速精細(xì)積分方法精細(xì)積分方法基于兩個(gè)要點(diǎn)，一個(gè)是加法定理，另一個(gè)是增量計(jì)算和存儲(chǔ)。 對(duì)于給定矩陣h,它對(duì)應(yīng)的矩陣指數(shù)exp(h)有如下性質(zhì)，即exp(h) = expj2arexp(h92'，hzdi)如果tv足夠大，則矩陣it比較小，那么矩陣it的矩陣指數(shù)可用如下的階taylor 級(jí)數(shù)近似，即q'.即(12)<d = exp(h)-i + h +精細(xì)積分方法s7將h'的矩陣指數(shù)分為兩部分，中-i + r， r=h +，一(13)然后對(duì)

7、增量部分應(yīng)用加法定理，即r = r2+2r(14)通過(guò)執(zhí)行tv次方程(14),然后將r加上單位矩陣，則得到h對(duì)應(yīng)的矩陣指數(shù)。上面簡(jiǎn)要地介紹了精細(xì)積分方法，此方法編寫(xiě)程序，計(jì)算精度非常高。但是, 對(duì)于大規(guī)模系統(tǒng)，由于系統(tǒng)的自由度數(shù)目很大，計(jì)算矩陣指數(shù)將非常耗費(fèi)時(shí)間和 內(nèi)存。雖然，對(duì)于大規(guī)模動(dòng)力系統(tǒng)，其剛度、阻尼和質(zhì)量矩陣是稀疏矩陣，從而 h也是稀疏矩陣，但是直接通過(guò)以上的精細(xì)積分方法計(jì)算其矩陣指數(shù)，在計(jì)算過(guò)程，特別是方程(14)的加法定理的合并過(guò)程中，矩陣將逐漸變?yōu)闈M(mǎn)矩陣或非常稠 密的矩陣，很難利用矩陣的稀疏性質(zhì)，從而計(jì)算量很大。本文利用結(jié)構(gòu)動(dòng)力響應(yīng)的物理特點(diǎn)，從物理上說(shuō)明，大規(guī)模動(dòng)力系統(tǒng)對(duì)應(yīng)的

8、 矩陣指數(shù)理論上應(yīng)該是稀疏矩陣，這樣在計(jì)算過(guò)程中可大大節(jié)約計(jì)算量，從而給 出一種計(jì)算矩陣指數(shù)的高效算法，且可節(jié)省存儲(chǔ)要求。而原始精細(xì)積分方法之所 以得到非常稠密的矩陣是因?yàn)橛?jì)算誤差造成的。眾所周知的事實(shí)是能量在一維桿屮的傳播速度是有限值同理，在離散的結(jié)構(gòu)中，雖然其能量的傳播速度很難確切得到，但其動(dòng)力學(xué)響應(yīng)的能量傳遞 速度肯定是有限的，這將對(duì)矩陣指數(shù)的結(jié)構(gòu)產(chǎn)生很大影響。根據(jù)方程(6)，如果 外力為零，貝ij方程(6)可寫(xiě)為如下分塊形式，即y+l = t.y<- + ti2p(15p+l= t2iya- + t22pa-由上述方程可以得到矩陣指數(shù)元素的物理意義，g卩：tu的第f行第列元素表示

9、 第y個(gè)自由度上給定單位位移，而其他&由度位移為零且所有自由度動(dòng)量為零 時(shí)，經(jīng)過(guò)一個(gè)積分步長(zhǎng)77后，第/個(gè)自由度的位移響應(yīng)；t12的第/行第j列元素 表示第/個(gè)自由度上給定單位動(dòng)量，而艽他自由度動(dòng)量為零且所有自由度位移為 零時(shí)，經(jīng)過(guò)一個(gè)積分步長(zhǎng)7后，第/個(gè)自由度的位移響應(yīng)；t2i的第第y列元 素表示第./個(gè)自由度上給定單位位移，而其他自由度位移為零且所有自由度動(dòng)量 為零時(shí)，經(jīng)過(guò)一個(gè)積分步長(zhǎng)7?后，第/個(gè)自由度的動(dòng)量響應(yīng)；t22的第/行第y列元素表示第j個(gè)自由度上給定單位動(dòng)量，而其他自由度動(dòng)量為零且所有自由度位移為零吋，經(jīng)過(guò)一個(gè)積分步長(zhǎng)z/后，第/個(gè)自由度的動(dòng)量響應(yīng)。由于能量在結(jié)構(gòu)屮的傳

10、遞速度是有限的，假設(shè)某個(gè)自由度上有擾動(dòng)，在較小 的吋間內(nèi)，必然只能傳播到有限的自由度，而不可能傳播到所有自由度。根據(jù)上 面給出的丁的物理含義，則它們一-定是稀疏矩陣。這樣，既可以將矩陣指數(shù)作為稀疏矩陣計(jì)算和存儲(chǔ)，從而節(jié)約計(jì)算量和存儲(chǔ)空間。對(duì)于給定矩陣h和積分步長(zhǎng)/z,原始的精細(xì)積分方法按如k過(guò)程計(jì)算矩陣指 數(shù)：首先確定然后令it = h|,其次按照方程(13)計(jì)算r,然后執(zhí)行n次方程(14),最后將r與單位矩陣相加得到矩陣指數(shù)。若采用上述的計(jì)算過(guò)程，雖然矩陣h為稀疏矩陣，但是如果觀察根據(jù)方程(13)計(jì)算得到的矩陣r，可以發(fā)現(xiàn)它可能是一個(gè)滿(mǎn)矩陣或很不稀疏的矩陣，但仔細(xì)觀 察會(huì)發(fā)現(xiàn)，此矩陣有很多很

11、小的元素。另一方面，根據(jù)上面的分析，此時(shí)的矩陣 r相當(dāng)于很小的時(shí)間步長(zhǎng)|對(duì)應(yīng)的矩陣指數(shù)減去*個(gè)單位矩陣，因此按照上面分析的能量傳播的性質(zhì)，它一定是非常稀疏的矩陣。因此，此時(shí)r矩陣中非常小 的非零元素是計(jì)算屮的誤差，它們理論上應(yīng)該為零，應(yīng)該將它們判斷出來(lái)，并令 它們?yōu)榱?，從而將r轉(zhuǎn)換為稀疏矩陣存儲(chǔ)。具體過(guò)程如下：從上面分析給出的矩 陣指數(shù)的物理意義可知，r矩陣的四塊子矩陣的物理意義不同，因此我們將r分 為，、，四塊，假設(shè)為矩陣ru中絕對(duì)值最大的元素，并給定一個(gè)誤 差標(biāo)準(zhǔn)，如e=l(t25,則rn中的元素如果滿(mǎn)足其絕對(duì)值小于e' h，則認(rèn)為此元素應(yīng)該為零，根據(jù)此原則可將rh稀疏存儲(chǔ)。同樣，

12、可將ri2，r2i，r22稀疏化。以上過(guò)程定義為一個(gè)矩陣的稀疏化。同理，方程(14)給出的加法定理同樣有上述 類(lèi)似的現(xiàn)象，因此，若直接應(yīng)用方程(14),兒次合并后，矩陣將變?yōu)闈M(mǎn)陣或稠密 矩陣，同樣可進(jìn)行上述的矩陣稀疏化。對(duì)于大規(guī)模動(dòng)力系統(tǒng)，對(duì)于一個(gè)合理的時(shí)間步長(zhǎng)，某一處的擾動(dòng)經(jīng)過(guò)一個(gè)時(shí) 間步長(zhǎng)后，其影響只是局部化的，不會(huì)傳播到整個(gè)結(jié)構(gòu)，因此系統(tǒng)對(duì)應(yīng)的矩陣指 數(shù)一定是稀疏矩陣，則可將精細(xì)積分方法作如下改進(jìn)，即對(duì)于給定矩陣h和積分 步長(zhǎng)/7,我們可給出如下的快速精細(xì)積分算法(fpim)，來(lái)計(jì)算矩陣指數(shù)。1. 由于h是稀疏矩陣，將h按照稀疏矩陣存儲(chǔ)；2. 確定/v，y8,91;3. 計(jì)算 h、h|;

13、4. 計(jì)算 r = h' +;5. 將矩陣r稀疏化；6. 執(zhí)行如下語(yǔ)句，即for iter=l:nr=r*( r + 2*1);將r稀疏化；end7. 得到矩陣r后，將其與單位矩陣相加，即得到h和a對(duì)應(yīng)的矩陣指數(shù)。上述計(jì)算過(guò)程與原始精細(xì)積分方法相比，只是增加了矩陣r的稀疏化過(guò)程， 但是這樣的處理將極大地提高計(jì)算效率，具體比較請(qǐng)見(jiàn)數(shù)值算例。得到矩陣指數(shù) 后，還須考慮外力作用的影響。根據(jù)上面分析矩陣指數(shù)結(jié)構(gòu)的相同思想，可以對(duì) 外力部分采用完全相同的處理方法，由于基本思想完全相同，本文不作詳細(xì)介紹。4數(shù)值算例算例1:考慮如圖1所示的彈簧質(zhì)量系統(tǒng)，若取 =2000,則系統(tǒng)包貪2000個(gè) 質(zhì)點(diǎn)，

14、本文隨機(jī)選擇每個(gè)質(zhì)點(diǎn)的質(zhì)量和彈簧的剛度，結(jié)果分別如圖2和3。此結(jié) 構(gòu)很容易寫(xiě)出其剛度矩陣k和質(zhì)量矩陣m，而阻尼矩陣取為c=0.05k。假設(shè) 每個(gè)質(zhì)點(diǎn)上都作用相同的外力cos(20,則系統(tǒng)的運(yùn)動(dòng)方程為mjjfe- c和 ky = f 爵i(2/)+ cos(2z)(16)其中f = 1, 1, l , if(17)分別采用木文方法(fpim)、原始精細(xì)積分方法(pim)、4階r-k方法和 newmark方法積分此問(wèn)題，稅分區(qū)間為0- 1000(s)o本文方法和原始精細(xì)稅分方法的積分步長(zhǎng)為l(s), 4階r-k方法采用matlab的ode45函數(shù)計(jì)算，并將絕對(duì) 誤差和相對(duì)誤差均設(shè)置為10_ 13,

15、 newmark方法的積分步長(zhǎng)為0.001(s)。以r-k方法為參考解，分別計(jì)算木文方法、原始精細(xì)積分方法和newmark方法的相對(duì) 誤差。位移和動(dòng)量的相對(duì)誤差分別定義為edl|y- yrklkiirkl|p- prk|2tpt(18)其屮yrk和分別表示r-k方法給出的位移和動(dòng)量，而y和p可取本文方法、原始精細(xì)積分方法或newmark方法積分得到的位移和動(dòng)量。本文方法積分到1000(s)時(shí)，各個(gè)質(zhì)點(diǎn)的位移和動(dòng)暈如圖4所示。newmark方法、本文方法和原始精細(xì)積分方法的相對(duì)誤差如圖5中。圖5表明，本文方法 和原始精細(xì)積分方法的精度都非常高，而本文方法和原始積分方法給出的結(jié)果非 常接近，并沒(méi)有

16、損失計(jì)算精度，這說(shuō)明了本文方法的正確性。圖5同時(shí)表明 newmark方法采用0.001(s)步長(zhǎng)積分，精度也沒(méi)右達(dá)到本文方法的精度。四種方法的計(jì)算時(shí)間如表1所示，可以看到，本文方法的計(jì)算效率要大大優(yōu) 于原始的精細(xì)積分方法、matlab的ode45和newmark方法。對(duì)于matlab的 ode45和newmark方法要達(dá)到較高的計(jì)算精度，必須采用小的多的積分步長(zhǎng)， 因此效率降低，newmark方法要達(dá)到更高的精度，還必須減小步長(zhǎng)，從而效率 還耍降低。按照前文的分析，由于矩陣指數(shù)中存在大量的零元素，原始精細(xì)積分 方法效率較低的原因在于浪費(fèi)了大量零元素的乘法運(yùn)算。圖1彈簧質(zhì)量系統(tǒng)圖2各質(zhì)點(diǎn)的質(zhì)量

17、luauomidst;卜 -1ur. fi*，k-l»-?l:7 6 5 4 3 2 ? z?c.s»601801 # 10011201140116011801*2001node圖3各彈簧的剛度b)5 0 5e3bcoe60180110011201140116011801node60180110011201140116011801node圖4本文方法給出的(=1000(s)時(shí)各質(zhì)點(diǎn)的位移和動(dòng)量iou10叫neumarkio4110叫neumark10義圖5算例1相對(duì)誤差比較表1算例1計(jì)算效率比較fpimode45pimnewmarkcpu time(s) 18.31162.

18、7412.5440.4算例2:考慮如圖6所示的平而應(yīng)力的動(dòng)力學(xué)問(wèn)題，結(jié)構(gòu)由兩種材料組成，兩種 材料的密度和泊松比相同，而楊氏模量不同，它們分別為r = 8? 103(kg/m3), m 0.2, £, = 3? 1011 (n/m), £22? 10n(n/m) (19)有限元網(wǎng)格和節(jié)點(diǎn)編號(hào)如圖7所示，采用三節(jié)點(diǎn)三角形單元，質(zhì)量矩陣采用集中 質(zhì)量矩陣，共有3200個(gè)單元，1681個(gè)節(jié)點(diǎn)，施加邊界條件后，共有3280個(gè)自 由度。初始位移為零，各自由度上具有初始動(dòng)量1,無(wú)外力作用。分別采用木文方法(fpim)、原始精細(xì)積分方法(pim)、4階r-k方法和 newmark方法積、

19、分此問(wèn)題，積分區(qū)間為0- 4? 10' 3(s)o本文方法和原始精細(xì)釈 分方法的積分步長(zhǎng)為10_6(s)，4階r-k方法采用matlab的ode45函數(shù)計(jì)算，并將絕對(duì)誤差和相對(duì)誤差均設(shè)置為10_ 13, newmark方法的積分步長(zhǎng)為10_ 8(s)。以r-k方法為參考解，分別計(jì)算木文方法、原始精細(xì)積分方法和newmark方法 的相對(duì)誤差。位移和動(dòng)量的相對(duì)誤差的定義與算例1相同。本文方法積分到4' i(y30o時(shí)，各個(gè)質(zhì)點(diǎn)x方向的位移和動(dòng)量如圖8所示，而方向的位移和動(dòng)量如圖9所示。newmark方法、本文方法和原始精細(xì)積分方法的相對(duì)誤差如圖10。圖10表明，本文方法和原始精細(xì)積

20、分方法的精度都很 高，并且本文方法和原始積分方法給出的結(jié)果非常接近，這再次說(shuō)明了本文方法 的正確性。圖10同時(shí)表明newmark方法采用10_ 8(s)步長(zhǎng)積分，相對(duì)誤差的精 度僅達(dá)到約10"3 (s)量級(jí)。四種方法的計(jì)算時(shí)間如表2所示，可以看到，本文方法的計(jì)算效率要優(yōu)于 matlab的ode45和newmark方法，比原始精細(xì)積分方法更是快了約220倍，與 原始精細(xì)積分方法相比效率得到了巨人的提高。lme'，r，me2, r, m0.5m 0.5mml <141'.*；*mvh2li4 - a 4« -m< 略“ ktjbtib4b461*2h

21、mmobwkmh5?ku嶙門(mén) i »u蟮<> 圖6兩種材料組成的平囬應(yīng)力m題hb w>«ui4 *>（«««>吻 ho參 ufl uhm«4»1»m4k，叫im1i >賽 >肩*|»«7*«魯 1km 叫< * *«瓤4*mh*mi 鴒 m“詹<«）和鴇 m k嗱 m磷 m 5b *m» * *钃 * » *«*|ha> »wmi m« umma似籲

22、1;縐吻 m».4 鷀林霣 mb m* mmr >» w» m靜籲如和礞>« w»| mv ml h k林*mfl m v4»k«« ?*<鉍摩 “h*m<<«hb5b*<?* mu him3i耗篇”魯？'mrmm 勿h 科鴒卿<» *鑛 ?<>，> k碥肇mt hu u>笮*蝙”鯫 uh吻 *«詹，wi4-m，mmx w414-w m«4e« * mp，廖m« >, ;ai

23、w如” , *#< 一 螬 ihi3«i|i we*-iu纗鴦*j»*h»-钃 ?蠊？wk9ir 蟾r k m霸mo*?瓣> 鑲 m > m« hh<om na m*m籌圖7有限元網(wǎng)格和節(jié)點(diǎn)編號(hào)圖84? l（rs）時(shí)刻各節(jié)點(diǎn)x方向的位移和動(dòng)量 0123隹 s!»oj;pxjo ilemuolu圖94? l（rs）時(shí)刻各節(jié)點(diǎn)y方向的位移和動(dòng)量圖10算例2相對(duì)誤差比較表2算例2計(jì)算效率比較fpimode45pimnewmarkcpu time(s) 105.51070.923252.15319.2結(jié)論本文提出了一種針對(duì)大規(guī)模

24、動(dòng)力系統(tǒng)的改進(jìn)的快速精細(xì)積分方法(fpim),利 用大規(guī)模動(dòng)力系統(tǒng)矩陣的稀疏性和動(dòng)力問(wèn)題能量傳遞有限的物理特性，表明對(duì)于 合理的積分步長(zhǎng)，矩陣指數(shù)具有稀疏結(jié)構(gòu)，并基于此給出一種計(jì)算大規(guī)模動(dòng)力系 統(tǒng)動(dòng)力響應(yīng)的高效率方法。數(shù)值算例表明，本文方法僅需在原始精細(xì)積分方法基 礎(chǔ)上進(jìn)行簡(jiǎn)單修改，即可將原始精細(xì)積分方法的效率大幅度提高，并且效率也比 常用的4階r-k方法和newmark方法高。參考文獻(xiàn)1. n. m. newmark. a method of computation for structural dynamics. ascejournal of engineering mechanics，8

25、5(3): 67-94, 1959.2. hairer e， lubich c， wanner g. geometric numerical integration:structure-preserving algorithm for ordinaiy differential equations (secondedition). new york: springer, 2006.3. hairer e, n0rsett s p，wanner g. solving ordinary differential equationsi-no ns tiff problems (second edit

26、ion). berlin: springer, 1993.4. hairer e，wanner g. solving ordinary differential equations ll-stiff anddifferential- algebraic problems (second edition). berlin: springer，1996.5. 鐘萬(wàn)勰。結(jié)構(gòu)動(dòng)力方程的精細(xì)時(shí)程積分法。大連理工大學(xué)學(xué)報(bào)，34(2):131-136， 1994。6. zhong wx, williams fw. a precise time step integration method. proceed

27、ings of the institution of mechanical engineers. part c. journal of mechanical engineering science 1994; 208(6): 427-430.7. zhong wx. on precise integration method. journal of computational and applied mathematics 2004; 163(1): 59-78.8. zhang hw, chen bs，gu yx. an adaptive algorithm of precise integration for transient analysis. acta mechannica solida sinica 2001; 14(3): 215-224.9. 譚述君，吳志剛，鐘萬(wàn)勰。矩陣指數(shù)精細(xì)積分方法中參數(shù)的自適應(yīng)選擇。力學(xué)學(xué)報(bào)，41(6)：961-966, 2009。10. 孫雁。奇異系統(tǒng)矩陣的精細(xì)積