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

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

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

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

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

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

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

1 改進(jìn)精細(xì)積分單步法

非線性動(dòng)力學(xué)方程的一般可表達(dá)成如下形式

(1)

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

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

(2)

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

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

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

(3)

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

i=1,2,3,…,s

構(gòu)造如下的隱式單步塊方法

(4)

將式(4)代入標(biāo)準(zhǔn)測試方程

(5)

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

(6)

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

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

(7)

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

(8)

(9)

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

(10)

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

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

(11)

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

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

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

(12)

當(dāng)s=3時(shí)，按照上述的構(gòu)造方法可以獲得相應(yīng)的B、d的表達(dá)式

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

(13)

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

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

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

(14)

式中：ci=i/4,i=0,1,3,4；A0=A4=1/18；A1=A3=4/9。Rk[f]是新四點(diǎn)積分公式的截?cái)嗾`差?？梢?，新四點(diǎn)積分公式的系數(shù)具有如下特征：

不難驗(yàn)證當(dāng)k=0,1,2,3時(shí)，Rk[f]≡0，當(dāng)k=4時(shí)，R4[f]=1/480。因此，式(13)具有3階代數(shù)精度，且在任意積分區(qū)間[a,b]上截?cái)嗾`差為

(15)

為了對(duì)比，給出同階三點(diǎn)和四點(diǎn)Newton-Cotes公式的截?cái)嗾`差分別為

(16)

(17)

由此可見，新四點(diǎn)積分公式具有更小的截?cái)嗾`差系數(shù)絕對(duì)值，因而數(shù)值積分的精度要優(yōu)于同階的Newton-Cotes公式。

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

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

(18)

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

2 算例分析

算例1 二自由度線性動(dòng)力學(xué)方程

(19)

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

(20)

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

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

圖1 位移的絕對(duì)誤差曲線

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

(21)

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

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

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

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

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

(22)

(23)

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

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

表3 文獻(xiàn)[5]中單步法與本文單步法的計(jì)算時(shí)間比較

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

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

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

算例4 考慮平方非線性的二自由度動(dòng)力方程為

(24)

(25)

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

表4 非線性二自由度動(dòng)力學(xué)方程x、y計(jì)算結(jié)果對(duì)比

通過表4可知，本文方法在步長較小時(shí)有著與細(xì)分解相當(dāng)?shù)挠?jì)算精度，DQ4的計(jì)算精度與BM3相當(dāng)，繼續(xù)將仿真步長擴(kuò)大五倍，BM3仍能取得與細(xì)分解較為接近的結(jié)果。

3 結(jié) 論

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

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

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

(3) 本文提出基于單步塊方法的改進(jìn)精細(xì)積分單步法比預(yù)估校正-辛?xí)r間子域法、文獻(xiàn)[5]中的單步法具有更高的計(jì)算精度和效率，比文獻(xiàn)[6]中的單步法計(jì)算精度要高，計(jì)算效率理論上兩者相當(dāng)；

(4) 結(jié)合文獻(xiàn)[6]中方法和改進(jìn)精細(xì)積分單步法可以構(gòu)成自適應(yīng)顯式精細(xì)積分單步法，該方法可以應(yīng)用于結(jié)構(gòu)地震碰撞反應(yīng)分析[17-18]。

附錄A

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

(A.1)

Hm為m次多項(xiàng)式構(gòu)成的集合，則R(m,n)={R∶R=P(x)/Q(x);P(x)∈Hm,Q(x)∈Hn}，分子為m次多項(xiàng)式，分母為n次多項(xiàng)式(除去恒為零的元素)的有理分式的集合。Padé逼近法是從冪級(jí)數(shù)出發(fā)獲得有理逼近的一種十分簡潔而且非常有效的方法，其基本思想是對(duì)于一個(gè)給定的形式冪級(jí)數(shù)，構(gòu)造一個(gè)有理函數(shù)，使該有理函數(shù)的Taylor展開有盡可能多的項(xiàng)和原來的冪級(jí)數(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)所得的近似式。