改進后的動力學方程精細積分法

劉素娟，卿光輝

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

改進后的動力學方程精細積分法

劉素娟，卿光輝

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

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

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

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

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

1 結構動力方程的增維分塊精細積分法

1.1 指數矩陣的精細積分法[2]

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

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

其中，τ為步長。

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

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

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

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

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

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

為了計算T，有:

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

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

當該語句循環結束后，再執行：

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

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

一般n維線性定常結構的動力方程可以表示為：

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

η 為時間步長，在區間[tk，tk+1]內，式（11）可視為常系數微分方程：

再利用迭代求解：

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

2 改進后的動力學方程精細積分方法

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

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

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

將上式根據（7）式的思想可得出:

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

因此，所求線性定常結構的動力學方程的解的迭代表達式為:

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

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

3 實例分析

取自文獻[7]，求解的結構動力學方程為：

初始值為：

解析解為：

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

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

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

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

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

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

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

表2 相同步長下，不同步數的時間比較

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

4 結束語

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

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

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

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

[4]陳飚松，顧元憲.瞬態熱傳導方程的子結構精細積分方法[J].應用力學學報，2001，18(1):14-19.

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

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

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

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

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

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

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

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

book=7,ebook=138