非線性動力方程的一種改進精細積分單步方法

劉冬兵， 王 永， 李博文， 奕仲飛， 張 磊， 黎 慧

(1.攀枝花學院 數(shù)學與計算機學院，四川 攀枝花 617000； 2.國網上海市電力公司 特高壓換流站分公司,上海 201413；3.國網四川省電力公司內江供電公司，四川 內江 641000； 4.三峽大學 電氣與新能源學院，湖北 宜昌 443002)

鐘萬勰[1]提出精細積分法為非線性動力方程的時域分析提供一種單步顯式高精度數(shù)值算法。此外，研究人員基于精細積分法衍生出了針對Duhamel積分項不同的精細積分法衍生格式，其中具有代表性的有高斯精細法[2]、精細庫塔法[3]、精細積分多步法[4]、精細積分單步法[5-6]、高精度直接積分法[7]等等。

高斯精細積分法的計算精度雖然較高，但是高斯積分點不包括區(qū)間端點，對于弱形式求積元分析不再適用。精細庫塔法的精度有限，長時間的數(shù)值積分無法保持無阻尼振動系統(tǒng)的幅值不變，具有較大的數(shù)值累積誤差；精細積分多步法需要進行一次預估-校正，計算的本質是復合積分，缺點是數(shù)值積分系數(shù)中存在負數(shù)，對于初值比較敏感。對于精細積分單步法，王海波等采用梯形積分公式和Romberg算法逼近Duhamel積分項，考慮到梯形積分公式的代數(shù)精度且待求變量的預估值(二階Runge-Kutta法)精度較低，該算法整體精度還有待提高。王永等提出的精細積分單步法結合了精細積分法和微分求積法的各自特點，計算精度相對要高，能夠用來分析結構動力方程的時程響應。Li等的方法涉及矩陣求逆不在本文所討論的范圍。文獻[8]結合了廣義精細積分法和預估-校正法，整體計算精度和效率較王海波等要高。文獻[9]結合泰勒級數(shù)展開式和廣義精細積分法推導出了特定類型載荷時線性動力方程的逐步積分公式，具有任意階精度。

從數(shù)值積分的角度來看，文獻[10]分析了微分求積法直接用于分析結構動力學方程的代數(shù)精度階數(shù)，指出了非均勻網格要明顯優(yōu)于均勻網格。但是，從Duhamel積分項來看，采用均勻網格可以減少系數(shù)矩陣求解數(shù)量。王永等采用單步s級時域微分求積法的最后一個方程對Duhamel積分項進行數(shù)值積分，當s=4時所給出的數(shù)值積分格式中的加權系數(shù)存在負數(shù)，這將給數(shù)值積分的穩(wěn)定性帶來了不穩(wěn)定的風險，可能會放大數(shù)值積分的誤差并造成最終的計算結果失真。為此，文中基于單步塊方法提出了改進精細積分單步法。

單步塊方法(one-step block methods)是20世紀70年代由Shampine和Watts提出的一種求解常微分方程初值問題的數(shù)值方法[11]。具體地，對于單步s級塊方法將單個步長區(qū)間內的s個內點組成一個塊向量，從而由一個初始時刻的待求變量初值計算得到s個內點處的待求變量近似值。值得注意的是，單步塊方法比較類似于時域微分求積法和隱式Runge-Kutta方法的，但是二者之間的區(qū)別在于單步塊方法在計算中使用了初始時刻的函數(shù)值。

本文在王永等的基礎上，將單步塊方法與精細積分法相結合，采用s級的單步塊方法的第s個方程對Duhamel積分項進行數(shù)值積分，從而形成了一種改進精細積分單步方法。通過對線性動力學方程、非線性單擺和Van der Pol振子方程的數(shù)值仿真，并與預估校正-辛時間子域法、現(xiàn)有單步法以及隱式積分算法進行數(shù)值結果比較，結果表明文中所提出算法能夠對剛性較低或較高的非線性動力方程進行時程分析。

1 改進精細積分單步法

非線性動力學方程的一般可表達成如下形式

(1)

式中：H為n維常系數(shù)矩陣；r為非線性廣義外力項。

在任意積分區(qū)間[tk,tk+1]內，非線性動力方程(1)的解可一般表達為如下同解積分方程

(2)

式中：vk+1、vk分別表示待求向量在tk+1、tk時刻的值，且積分區(qū)間長度為Δt=tk+1-tk。等式(2)右邊第一項中的指數(shù)矩陣T=eHΔt可采用精細積分得到，而Duhamel積分項可以采用數(shù)值積分公式例如梯形公式、Simpson公式等等近似計算。本文采用s級的單步塊方法的第s個方程計算Duhamel積分項。

1.1 隱式單步塊方法構造數(shù)值積分公式

為了便于表述，考慮如下的一階常微分方程的初值問題

(3)

令c=(t-tk)/h,t∈(tk,tk+1)，從而將單步區(qū)間[tk,tk+1]正則化為c∈[0,1]，則式(3)可以改寫為

i=1,2,3,…,s

構造如下的隱式單步塊方法

(4)

將式(4)代入標準測試方程

(5)

式中，Re(λ)<0，令z=λh，則式(4)的穩(wěn)定性函數(shù)為

(6)

式中:es=[0 … 0 1]T∈Rs×1;Is為s維的單位矩陣。采用Padé逼近對穩(wěn)定性函數(shù)R(z)進行函數(shù)逼近，可以構造出A-穩(wěn)定、L-穩(wěn)定的計算格式。首先，給出如下的定理[12]。

定理1對于(k,j)-Padé逼近，若k≤j≤k+2，則此Padé逼近是A-穩(wěn)定的；若k

(7)

式中，φi是待定系數(shù)，將式(7)代入到式(4)中可得

(8)

(9)

對于(k,s)-Padé逼近(k

(10)

因此，式(9)有唯一的非零解。求得系數(shù)矩陣B和d后代入式(4)便可得到基于不同(k,s)-Padé逼近的計算格式。下面分別采用均勻網格和CGL網格 (Chebyshev-Gauss-Lobatto網格)構造不同的計算格式進行，并對其第s個方程的數(shù)值積分的代數(shù)精度和截斷誤差進行分析。

考慮正則區(qū)間[0,1]上的均勻網格點分布

(11)

當s=3時，基于(2,3)-Padé逼近建立等式關系，可以求解得到B、d的表達式

同理，當?shù)胹=4時，基于(3,4)-Padé逼近構造計算格式的B、d系數(shù)矩陣表達式

考慮CGL網格在正則區(qū)間[0，1]上的網格點分布如下

(12)

當s=3時，按照上述的構造方法可以獲得相應的B、d的表達式

顯然，采用均勻網格的隱式單步塊方法的第s個方程可以用來作為數(shù)值積分公式，可以證明此時該數(shù)值積分公式就是著名的Newton-Cotes公式[13]。考慮采用CGL網格時得到新四點積分公式

(13)

式中：fk+i/4=f(tk+i/4,y(tk+i/4))，i=0,1,3,4。

1.2 新四點積分公式的數(shù)值精度及穩(wěn)定性

不失一般性，考慮一般多項式y(tǒng)′=xk在標準區(qū)間[0,1]上的積分值，并采用式(13)計算其近似值，則

(14)

式中：ci=i/4,i=0,1,3,4；A0=A4=1/18；A1=A3=4/9。Rk[f]是新四點積分公式的截斷誤差。可見，新四點積分公式的系數(shù)具有如下特征：

不難驗證當k=0,1,2,3時，Rk[f]≡0，當k=4時，R4[f]=1/480。因此，式(13)具有3階代數(shù)精度，且在任意積分區(qū)間[a,b]上截斷誤差為

(15)

為了對比，給出同階三點和四點Newton-Cotes公式的截斷誤差分別為

(16)

(17)

由此可見，新四點積分公式具有更小的截斷誤差系數(shù)絕對值，因而數(shù)值積分的精度要優(yōu)于同階的Newton-Cotes公式。

至此，利用式(13)完全可以采用王永等的思路實現(xiàn)式(2)中對Duhamel積分項的數(shù)值積分。

采用4階顯式龍格-庫塔法計算預估值yk+i/4(i=1,3,4)如下

(18)

針對文獻[6]不足，文中所提出的改進精細積分單步法，在計算量上改進精細積分單步法也只需要進行一次指數(shù)矩陣的精細積分即可，其余指數(shù)矩陣可通過矩陣乘法獲得。

2 算例分析

算例1 二自由度線性動力學方程

(19)

將其整理為式(1)的形式，則有：

(20)

式中，初始值x(0)=[0,0,0,0]T。

不同算法的仿真步長都取Δt=0.2 s，積分時程為10 s，預估公式都采用經典四階Runge-Kutta法，分別采用本文方法(簡記為BM3)、精細積分-微分求積法(簡記為DQ4)和精細科茨法(精細積分法加辛普森公式簡記為NC2，精細積分法加3階代數(shù)精度的科茨公式簡記為NC3)求解結果的相對于解析解的絕對誤差的絕對值R(x1)曲線如圖1所示。

圖1 位移的絕對誤差曲線

從圖1中可知，在代數(shù)精度相同的前提下，本文方法由于具有更小的截斷誤差系數(shù)絕對值，因而數(shù)值積分的精度更高。

(21)

式中：初值v1(0)=1.047 2；v2(0)=0。

以橢圓積分得到該問題的解析解為基準，分別采用本文方法(BM3)和精細積分-微分求積法(簡記為DQ4)計算幅角v1，計算步長取Δt=0.01 s，相應的計算結果列入表1中。

表1 幅角v1的數(shù)值結果對比(Δt=0.01 s)

從表1可知，本文方法的計算精度要明顯高于精細積分-微分求積法，且與橢圓積分解完全吻合，這再次驗證了本文方法在計算精度上的優(yōu)勢。

算例3 考慮Van der Pol振子方程為

(22)

(23)

表2 Van der Pol方程中x的計算結果對比(Δt=0.1 s)

從表2可知，本文方法的計算結果可以與細分解保持小數(shù)點后六位保持一致，計算精度要明顯高于同階的精細積分-微分求積法和預估校正-辛時間子域法[14]。表3是此時本文算法和文獻[5]中算法的計算時間對比，從中可知本文算法的計算效率較高，這種優(yōu)勢隨著計算步數(shù)的增加表現(xiàn)得更為突出。

表3 文獻[5]中單步法與本文單步法的計算時間比較

圖2 Van der Pol振子位移軌跡(BM3)

圖3 Van der Pol振子位移軌跡(TR)

從圖2～3可知，對于非線性剛性Van der Pol振子方程，采用本文方法也可以獲得比較精確的數(shù)值解。這說明雖然本文方法嚴格上屬于顯式積分算法，但是相對于更適合求解剛性微分方程的隱式梯形算法依然具有較大的優(yōu)勢。

算例4 考慮平方非線性的二自由度動力方程為

(24)

(25)

取待求變量的初值v1(0)=v2(0)=0.1，v3(0)=v4(0)=0。本文方法的步長分別在Δt=0.1 s和Δt=0.5 s時，細分解和本文方法以及DQ4的數(shù)值結果比較如表4所示。

表4 非線性二自由度動力學方程x、y計算結果對比

通過表4可知，本文方法在步長較小時有著與細分解相當?shù)挠嬎憔龋珼Q4的計算精度與BM3相當，繼續(xù)將仿真步長擴大五倍，BM3仍能取得與細分解較為接近的結果。

3 結 論

本文是對文獻[6]中思路的延續(xù)和改進，基于單步塊方法和精細積分法得出了如下的結論：

(1) 采用均勻網格的單步塊方法構造的數(shù)值積分公式即是Newton-Cotes公式；

(2) 采用CGL網格的單步塊方法構造的四點數(shù)值積分公式比同階代數(shù)精度的Newton-Cotes公式的數(shù)值積分精度要高；

(3) 本文提出基于單步塊方法的改進精細積分單步法比預估校正-辛時間子域法、文獻[5]中的單步法具有更高的計算精度和效率，比文獻[6]中的單步法計算精度要高，計算效率理論上兩者相當；

(4) 結合文獻[6]中方法和改進精細積分單步法可以構成自適應顯式精細積分單步法，該方法可以應用于結構地震碰撞反應分析[17-18]。

附錄A

Padé逼近是一種有理數(shù)逼近，它克服了用多項式逼近大撓度函數(shù)效果不理想，而用冪函數(shù)(如Taylor級數(shù))逼近函數(shù)收斂性太差等缺點。記

(A.1)

Hm為m次多項式構成的集合，則R(m,n)={R∶R=P(x)/Q(x);P(x)∈Hm,Q(x)∈Hn}，分子為m次多項式，分母為n次多項式(除去恒為零的元素)的有理分式的集合。Padé逼近法是從冪級數(shù)出發(fā)獲得有理逼近的一種十分簡潔而且非常有效的方法，其基本思想是對于一個給定的形式冪級數(shù)，構造一個有理函數(shù)，使該有理函數(shù)的Taylor展開有盡可能多的項和原來的冪級數(shù)相吻合。其具體定義如下

定義：如果有理函數(shù)R=P(x)/Q(x)滿足f(x)-P(x)/Q(x)=Ο(xm+n+1),Q(0)=1，則稱P(x)/Q(x)為f(x)在R(m,n)中的Padé逼近，記為(m,n)f。

依據(jù)Padé逼近的定義，f(x)的(m,n)-Padé逼近也可以理解為由方程

qm(x)f(x)-pn(x)=Ο(xm+n+1)

(A.2)

解出f(x)所得的近似式。