劉素娟，卿光輝

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

改進后的動力學(xué)方程精細積分法

劉素娟，卿光輝

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

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

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

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

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

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

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

按常微分方程理論，對于一個齊次方程：

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

其中，τ為步長。

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

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

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

文獻[5]指出指數(shù)截斷階數(shù)L=4，迭代次數(shù)N=20時精細積分法中指數(shù)矩陣計算可以認為很接近其精確解答。此時指數(shù)矩陣T與單位矩陣In相差不遠，因此可寫為：

其中，Ta陣是一個小量矩陣。

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

為了計算T，有:

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

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

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

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

1.2 增維分塊精細積分法[7]

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

對（9）式進行增維處理可得：

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

再利用迭代求解：

其中，exp（Akη）采用精細積分計算。

2 改進后的動力學(xué)方程精細積分方法

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

簡單檢驗可知矩陣Rk必為可逆的，其中

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

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

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

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

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

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

3 實例分析

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

初始值為：

解析解為：

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

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

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

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

圖1 新算法的x1與其準確解比較

圖2 新算法的x2與其準確解比較

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

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

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

4 結(jié)束語

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

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

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

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

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

[5]汪夢甫，區(qū)達光.精細積分方法的評估與改進[J].計算力學(xué)學(xué)報，2004，21(6):728-733.

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

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

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

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

[10]張慶云，滕圣剛.一種提高增維精細積分法計算精度的方法[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

中國民航大學(xué)校級科研基金項目(編號：2012kye07)

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

book=7,ebook=138