改進(jìn)后的動(dòng)力學(xué)方程精細(xì)積分法

劉素娟，卿光輝

（中國(guó)民航大學(xué)航空工程學(xué)院，天津 300300）

改進(jìn)后的動(dòng)力學(xué)方程精細(xì)積分法

劉素娟，卿光輝

（中國(guó)民航大學(xué)航空工程學(xué)院，天津 300300）

對(duì)于結(jié)構(gòu)動(dòng)力學(xué)方程，在分塊增維精細(xì)積分法的基礎(chǔ)上，使各個(gè)子塊產(chǎn)生共同項(xiàng)，即子塊相似化。這樣精細(xì)積分每次循環(huán)的結(jié)果，多個(gè)子塊可以同時(shí)利用，同時(shí)也減少了計(jì)算量和存儲(chǔ)量。在非齊次項(xiàng)的處理問題上，取每個(gè)步長(zhǎng)的中間值，簡(jiǎn)化了計(jì)算，也保留了精度。數(shù)值例題顯示研究方法的高效性和可行性。

動(dòng)力學(xué)方程；增維分塊；精細(xì)積分法；子塊矩陣；相似化

精細(xì)積分問世以來，結(jié)構(gòu)動(dòng)力學(xué)方程求解的過程大為簡(jiǎn)化。Moler and Van Loan[1]提到了十九種計(jì)算方法，并就多種方法的穩(wěn)定性、適應(yīng)性和計(jì)算精度等問題作了較詳細(xì)的分析和討論。鐘萬勰院士等人提出的精細(xì)積分法[2]，在應(yīng)用于求解線性定常結(jié)構(gòu)動(dòng)力方程時(shí)，能夠在數(shù)值上得到逼近于機(jī)器精度的結(jié)果，在相關(guān)領(lǐng)域[3～6]得到了廣泛應(yīng)用。而對(duì)于非齊次的結(jié)構(gòu)動(dòng)力方程，顧元憲、陳飚松、張洪武等人提出的增維精細(xì)積分法[7]避免了矩陣求逆，同時(shí)提高了數(shù)值計(jì)算的穩(wěn)定性。但是增維精細(xì)積分方法一般要求每一個(gè)時(shí)間步都要重新通過20次矩陣加、乘迭代來計(jì)算T矩陣，所需的計(jì)算量和計(jì)算時(shí)間將很大，特別是在處理大型問題時(shí)。張繼鋒，鄧子辰提出的結(jié)構(gòu)動(dòng)力方程的增維分塊精細(xì)積分法[9]，在增維精細(xì)積分法的基礎(chǔ)上，對(duì)矩陣進(jìn)行分塊計(jì)算，考慮非齊次項(xiàng)的特點(diǎn)，減小了矩陣的維數(shù)，實(shí)現(xiàn)簡(jiǎn)化計(jì)算，但是在計(jì)算速度和非齊次項(xiàng)處理上都不是很理想。張慶云、滕圣剛提出的一種提高增維精細(xì)積分法計(jì)算精度的方法[10]，在每一個(gè)時(shí)間步長(zhǎng)內(nèi)，仍然將非齊次項(xiàng)當(dāng)成常數(shù)，但是該常數(shù)的值取為該時(shí)間段內(nèi)不同時(shí)刻值的平均值，或者取為中間時(shí)刻的值，計(jì)算精度得到了一定的改善，但是整體計(jì)算量的減少并不理想。

本文在前人所做工作的基礎(chǔ)上，借鑒已有增維分塊精細(xì)積分法的思想，對(duì)不同子塊進(jìn)行化簡(jiǎn)處理，使各個(gè)子塊產(chǎn)生共同項(xiàng)，或者產(chǎn)生每次循環(huán)結(jié)果項(xiàng)，即子塊結(jié)構(gòu)相似化處理。這樣不僅可以使相同的項(xiàng)一次性計(jì)算，還能使精細(xì)積分每次循環(huán)的結(jié)果，多個(gè)子塊可以同時(shí)利用，同時(shí)也減少了計(jì)算量和存儲(chǔ)量，提高了計(jì)算速度。另外，在非齊次項(xiàng)常數(shù)化處理時(shí)，本文采用最新的張慶云、滕圣剛[10]的非齊次項(xiàng)取值為該時(shí)間段內(nèi)不同時(shí)刻值的平均值，或者中間時(shí)刻的值，使計(jì)算精度大為提高。因而本文在計(jì)算速度和計(jì)算精度上都有很大的優(yōu)勢(shì)，尤其是計(jì)算步數(shù)較大的問題時(shí)優(yōu)勢(shì)更加明顯。

1 結(jié)構(gòu)動(dòng)力方程的增維分塊精細(xì)積分法

1.1 指數(shù)矩陣的精細(xì)積分法[2]

按常微分方程理論，對(duì)于一個(gè)齊次方程：

H為定常矩陣，其通解可以表示為:

其中，τ為步長(zhǎng)。

精細(xì)積分法主要是將注意力放在增量上，而不是全量。假設(shè)A為n×n階矩陣，為了避免舍入誤差，利用指數(shù)函數(shù)的加法定理：

其中，m為任意正整數(shù)，且m=2N（例如N=20，則m=1 048 576)。

由于τ是一個(gè)小的時(shí)間區(qū)段，則η=τ/m是一個(gè)非常小的時(shí)間區(qū)段。因此，對(duì)于η，有

文獻(xiàn)[5]指出指數(shù)截?cái)嚯A數(shù)L=4，迭代次數(shù)N=20時(shí)精細(xì)積分法中指數(shù)矩陣計(jì)算可以認(rèn)為很接近其精確解答。此時(shí)指數(shù)矩陣T與單位矩陣In相差不遠(yuǎn)，因此可寫為：

其中，Ta陣是一個(gè)小量矩陣。

因?yàn)門a很小，當(dāng)它與單位陣In相加時(shí)，就會(huì)成為尾數(shù)，在計(jì)算機(jī)的舍入操作中會(huì)被舍去，影響計(jì)算精度[8]。所以，至關(guān)重要的一點(diǎn)是指數(shù)矩陣的存儲(chǔ)是式(3)中的 Ta，而不是 Ta+In。

為了計(jì)算T，有:

這種分解一直做下去，共N次。其次應(yīng)注意，對(duì)任意矩陣Tb和Tc有:

將Tb和Tc都看成為Ta，因此式（4）的 N次乘法在計(jì)算機(jī)中相當(dāng)于以下的語句：

當(dāng)該語句循環(huán)結(jié)束后，再執(zhí)行：

執(zhí)行N次乘法后，Ta已不再是很小的矩陣了，這個(gè)加法已沒有嚴(yán)重的舍入誤差。

1.2 增維分塊精細(xì)積分法[7]

一般n維線性定常結(jié)構(gòu)的動(dòng)力方程可以表示為：

對(duì)（9）式進(jìn)行增維處理可得：

η 為時(shí)間步長(zhǎng)，在區(qū)間[tk，tk+1]內(nèi)，式（11）可視為常系數(shù)微分方程：

再利用迭代求解：

其中，exp（Akη）采用精細(xì)積分計(jì)算。

2 改進(jìn)后的動(dòng)力學(xué)方程精細(xì)積分方法

對(duì)方程（12）進(jìn)行精細(xì)積分，其中N=20，對(duì)方程（17）進(jìn)行精細(xì)積分計(jì)算時(shí)，取N=20，每一個(gè)時(shí)間步都要重新通過20次矩陣加、乘迭代來計(jì)算矩陣,所需的計(jì)算量和計(jì)算時(shí)間將很大。借鑒文獻(xiàn)[5]的思想,對(duì)矩陣進(jìn)行分塊計(jì)算。

簡(jiǎn)單檢驗(yàn)可知矩陣Rk必為可逆的，其中

式中的兩個(gè)子塊都含有TaR，相似度很高，這是本文第一次成功的實(shí)現(xiàn)了子塊的相似處理。

將上式根據(jù)（7）式的思想可得出:

聯(lián)合（18）式和（19）式，可以看出矩陣T的子塊T11和 T12中都存在（I+TaR）2i，這使得兩個(gè)子塊存在很高的相似度，這是本文第二次成功的實(shí)現(xiàn)了子塊的相似化處理。尤其，當(dāng)T11采用精細(xì)積分求解，T11執(zhí)（7b）式時(shí)，T11從第1次到第N-1次，每一次環(huán)的結(jié)果，剛好可以被子塊T12中的（I+TaR）2i項(xiàng)完全利用。根據(jù)文獻(xiàn)[9]可知：增維后齊次方程的解的迭代表達(dá)式為：

因此，所求線性定常結(jié)構(gòu)的動(dòng)力學(xué)方程的解的迭代表達(dá)式為:

結(jié)合以上精細(xì)積分的全過程，只要確定初始值，就能很簡(jiǎn)單的計(jì)算出不同步長(zhǎng)和不同步數(shù)的方程的解。

非齊次項(xiàng)的處理采用文獻(xiàn)[10]的方法，將其值取為該時(shí)間段內(nèi)幾個(gè)不同時(shí)刻值的平均值，或者取為中間時(shí)刻的值，經(jīng)過計(jì)算分析發(fā)現(xiàn)，取中間時(shí)刻值滿足計(jì)算精度的同時(shí)，計(jì)算量較小。因此，本文取一個(gè)時(shí)間步長(zhǎng)[tk，tk+1]內(nèi)的中間時(shí)刻值，其表達(dá)式如下:

3 實(shí)例分析

取自文獻(xiàn)[7]，求解的結(jié)構(gòu)動(dòng)力學(xué)方程為：

初始值為：

解析解為：

選擇步數(shù)為1（第1步）的情況，步長(zhǎng)由0.01 s增加至0.50 s，每次增加0.01 s，用Mathematica語言編程得到x1，x2的解分別都是50組，并且將其與x1，x2的準(zhǔn)確解進(jìn)行對(duì)比，表1將列舉出50組數(shù)據(jù)中的8組數(shù)據(jù)。

表1 x1，x2新算法的解與準(zhǔn)確解其中8組數(shù)據(jù)的比較

比較分析這8組數(shù)據(jù)可知，本方方法跟準(zhǔn)確值相比，誤差很小，即使步長(zhǎng)增加至0.50，其結(jié)果也能精確至0.01位。再將這50組數(shù)據(jù)，用Mathematica軟件制成圖，橫軸為變化的步長(zhǎng)，縱軸為計(jì)算的結(jié)果，如圖1、圖2所示。

據(jù)圖1和圖2顯示并分析可知，在步長(zhǎng)由0.01增至0.50的整個(gè)過程，本文算法x1、x2的解與準(zhǔn)確解所代表的每組圖形幾乎完全重合，再次證明了本文方法的準(zhǔn)確性與高精度。

圖1 新算法的x1與其準(zhǔn)確解比較

圖2 新算法的x2與其準(zhǔn)確解比較

在此之前，計(jì)算結(jié)構(gòu)動(dòng)力學(xué)方程較快的方法是文[5]的增維分塊法，現(xiàn)在將本文方法與文[5]的方法的計(jì)算速度進(jìn)行對(duì)比。繼續(xù)分析上述算例，在步長(zhǎng)為0.02的條件下，當(dāng)步數(shù)分別為50 000、100 000、500 000、1 000 000、5 000 000、10 000 000，本文方法與文[5]增維分塊法計(jì)算時(shí)間對(duì)比，結(jié)果如表2所示。

表2 相同步長(zhǎng)下，不同步數(shù)的時(shí)間比較

由表2可以看出，這六組數(shù)據(jù)中本文方法的時(shí)間都比對(duì)對(duì)應(yīng)的文[5]的時(shí)間要短，當(dāng)計(jì)算的步數(shù)增加時(shí)，本文方法節(jié)省的時(shí)間就越多，因此計(jì)算的優(yōu)越性就越明顯，這對(duì)于解決大步數(shù)以及超大步數(shù)的結(jié)構(gòu)動(dòng)力學(xué)微分方程是有利的。

4 結(jié)束語

本文在增維分塊的基礎(chǔ)上，結(jié)合對(duì)子塊矩陣使各個(gè)子塊產(chǎn)生共同項(xiàng)，或者產(chǎn)生每次循環(huán)結(jié)果項(xiàng)，進(jìn)行結(jié)構(gòu)相似化處理的核心思想，提出了一種精確快速的，且程序上易實(shí)現(xiàn)解決結(jié)構(gòu)動(dòng)力學(xué)問題的快速精細(xì)積分方法。在步數(shù)很大時(shí)，不但能保證計(jì)算精度，而且計(jì)算速度更快，更能讓人滿意。

[1]MOLER C,LOAN C V.Nineteen dubiousways to compute the exponentialofamatrix[J].SIAM Rev，1978，20(4):801-836.

[2]鐘萬勰.結(jié)構(gòu)動(dòng)力方程的精細(xì)時(shí)程積分[J].大連理工大學(xué)學(xué)報(bào)，1994，34(2):131-136.

[3]顧元憲，陳飚松.瞬態(tài)熱傳導(dǎo)方程精細(xì)積分方法中對(duì)稱性的利用[J].力學(xué)與實(shí)踐，2000，22(5):19-22.

[4]陳飚松，顧元憲.瞬態(tài)熱傳導(dǎo)方程的子結(jié)構(gòu)精細(xì)積分方法[J].應(yīng)用力學(xué)學(xué)報(bào)，2001，18(1):14-19.

[5]汪夢(mèng)甫，區(qū)達(dá)光.精細(xì)積分方法的評(píng)估與改進(jìn)[J].計(jì)算力學(xué)學(xué)報(bào)，2004，21(6):728-733.

[6]候秀慧，鄧子辰，黃立新.橋梁結(jié)構(gòu)移動(dòng)荷載識(shí)別的辛精細(xì)積分算法[J].動(dòng)力學(xué)與控制學(xué)報(bào)，2008，6(1):66-72.

[7]顧元憲，陳飚松，張洪武.結(jié)構(gòu)動(dòng)力方程的增維精細(xì)積分法[J].力學(xué)學(xué)報(bào)，2000，32（4）：447-456.

[8]任 偉，杜鐵鈞.定常結(jié)構(gòu)動(dòng)力方程增維精細(xì)積分法求解的注記[J].杭州電子科技大學(xué)學(xué)報(bào)，2005，25(1):41-43.

[9]張繼鋒，鄧子辰.結(jié)構(gòu)動(dòng)力方程的增維分塊精細(xì)積分法[J].振動(dòng)與沖擊，2008，27(12):88-90.

[10]張慶云，滕圣剛.一種提高增維精細(xì)積分法計(jì)算精度的方法[J].科學(xué)技術(shù)與工程，2010，10(31):7627-7630.

The Kinetic Equation of Precise Integration Method Im proved

LIU Su-juan，QINGGuang-hui
(CivilAviation University ofChina College of Aeronautical Engineering，Tianjin 300300，China)

The structure equations,based on the block of increment dimensional precise integration method,so that each sub block into the common items,namely sub blocksimilarity.The precise integration of each cycle results,multiple sub blocks can beused at the same time,butalso reduce the computation and storage.To dealwith problems in nonhomogeneous term,intermediate values of each step,and simplifies the calculation,but also retains the accuracy.Numericalexamples show the efficiency and feasibility of researchmethods.

dynamic equation；incrementdimensionalprecise integrationmethod；block;blockmatrix；similarity

TU311.3

A

1672-545X（2014）04-0004-04

2014-01-07

中國(guó)民航大學(xué)校級(jí)科研基金項(xiàng)目(編號(hào)：2012kye07)

劉素娟(1987—)，女，湖北黃岡人，碩士，研究方向：復(fù)合材料力學(xué)。

book=7,ebook=138