非线性动力方程的一种改进精细积分单步方法

刘冬兵， 王 永， 李博文， 奕仲飞， 张 磊， 黎 慧

(1.攀枝花学院 数学与计算机学院，四川 攀枝花 617000； 2.国网上海市电力公司 特高压换流站分公司,上海 201413；3.国网四川省电力公司内江供电公司，四川 内江 641000； 4.三峡大学 电气与新能源学院，湖北 宜昌 443002)

钟万勰[1]提出精细积分法为非线性动力方程的时域分析提供一种单步显式高精度数值算法。此外，研究人员基于精细积分法衍生出了针对Duhamel积分项不同的精细积分法衍生格式，其中具有代表性的有高斯精细法[2]、精细库塔法[3]、精细积分多步法[4]、精细积分单步法[5-6]、高精度直接积分法[7]等等。

高斯精细积分法的计算精度虽然较高，但是高斯积分点不包括区间端点，对于弱形式求积元分析不再适用。精细库塔法的精度有限，长时间的数值积分无法保持无阻尼振动系统的幅值不变，具有较大的数值累积误差；精细积分多步法需要进行一次预估-校正，计算的本质是复合积分，缺点是数值积分系数中存在负数，对于初值比较敏感。对于精细积分单步法，王海波等采用梯形积分公式和Romberg算法逼近Duhamel积分项，考虑到梯形积分公式的代数精度且待求变量的预估值(二阶Runge-Kutta法)精度较低，该算法整体精度还有待提高。王永等提出的精细积分单步法结合了精细积分法和微分求积法的各自特点，计算精度相对要高，能够用来分析结构动力方程的时程响应。Li等的方法涉及矩阵求逆不在本文所讨论的范围。文献[8]结合了广义精细积分法和预估-校正法，整体计算精度和效率较王海波等要高。文献[9]结合泰勒级数展开式和广义精细积分法推导出了特定类型载荷时线性动力方程的逐步积分公式，具有任意阶精度。

从数值积分的角度来看，文献[10]分析了微分求积法直接用于分析结构动力学方程的代数精度阶数，指出了非均匀网格要明显优于均匀网格。但是，从Duhamel积分项来看，采用均匀网格可以减少系数矩阵求解数量。王永等采用单步s级时域微分求积法的最后一个方程对Duhamel积分项进行数值积分，当s=4时所给出的数值积分格式中的加权系数存在负数，这将给数值积分的稳定性带来了不稳定的风险，可能会放大数值积分的误差并造成最终的计算结果失真。为此，文中基于单步块方法提出了改进精细积分单步法。

单步块方法(one-step block methods)是20世纪70年代由Shampine和Watts提出的一种求解常微分方程初值问题的数值方法[11]。具体地，对于单步s级块方法将单个步长区间内的s个内点组成一个块向量，从而由一个初始时刻的待求变量初值计算得到s个内点处的待求变量近似值。值得注意的是，单步块方法比较类似于时域微分求积法和隐式Runge-Kutta方法的，但是二者之间的区别在于单步块方法在计算中使用了初始时刻的函数值。

本文在王永等的基础上，将单步块方法与精细积分法相结合，采用s级的单步块方法的第s个方程对Duhamel积分项进行数值积分，从而形成了一种改进精细积分单步方法。通过对线性动力学方程、非线性单摆和Van der Pol振子方程的数值仿真，并与预估校正-辛时间子域法、现有单步法以及隐式积分算法进行数值结果比较，结果表明文中所提出算法能够对刚性较低或较高的非线性动力方程进行时程分析。

1 改进精细积分单步法

非线性动力学方程的一般可表达成如下形式

(1)

式中：H为n维常系数矩阵；r为非线性广义外力项。

在任意积分区间[tk,tk+1]内，非线性动力方程(1)的解可一般表达为如下同解积分方程

(2)

式中：vk+1、vk分别表示待求向量在tk+1、tk时刻的值，且积分区间长度为Δt=tk+1-tk。等式(2)右边第一项中的指数矩阵T=eHΔt可采用精细积分得到，而Duhamel积分项可以采用数值积分公式例如梯形公式、Simpson公式等等近似计算。本文采用s级的单步块方法的第s个方程计算Duhamel积分项。

1.1 隐式单步块方法构造数值积分公式

为了便于表述，考虑如下的一阶常微分方程的初值问题

(3)

令c=(t-tk)/h,t∈(tk,tk+1)，从而将单步区间[tk,tk+1]正则化为c∈[0,1]，则式(3)可以改写为

i=1,2,3,…,s

构造如下的隐式单步块方法

(4)

将式(4)代入标准测试方程

(5)

式中，Re(λ)<0，令z=λh，则式(4)的稳定性函数为

(6)

式中:es=[0 … 0 1]T∈Rs×1;Is为s维的单位矩阵。采用Padé逼近对稳定性函数R(z)进行函数逼近，可以构造出A-稳定、L-稳定的计算格式。首先，给出如下的定理[12]。

定理1对于(k,j)-Padé逼近，若k≤j≤k+2，则此Padé逼近是A-稳定的；若k

(7)

式中，φi是待定系数，将式(7)代入到式(4)中可得

(8)

(9)

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

(10)

因此，式(9)有唯一的非零解。求得系数矩阵B和d后代入式(4)便可得到基于不同(k,s)-Padé逼近的计算格式。下面分别采用均匀网格和CGL网格 (Chebyshev-Gauss-Lobatto网格)构造不同的计算格式进行，并对其第s个方程的数值积分的代数精度和截断误差进行分析。

考虑正则区间[0,1]上的均匀网格点分布

(11)

当s=3时，基于(2,3)-Padé逼近建立等式关系，可以求解得到B、d的表达式

同理，当得s=4时，基于(3,4)-Padé逼近构造计算格式的B、d系数矩阵表达式

考虑CGL网格在正则区间[0，1]上的网格点分布如下

(12)

当s=3时，按照上述的构造方法可以获得相应的B、d的表达式

显然，采用均匀网格的隐式单步块方法的第s个方程可以用来作为数值积分公式，可以证明此时该数值积分公式就是著名的Newton-Cotes公式[13]。考虑采用CGL网格时得到新四点积分公式

(13)

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

1.2 新四点积分公式的数值精度及稳定性

不失一般性，考虑一般多项式y′=xk在标准区间[0,1]上的积分值，并采用式(13)计算其近似值，则

(14)

式中：ci=i/4,i=0,1,3,4；A0=A4=1/18；A1=A3=4/9。Rk[f]是新四点积分公式的截断误差。可见，新四点积分公式的系数具有如下特征：

不难验证当k=0,1,2,3时，Rk[f]≡0，当k=4时，R4[f]=1/480。因此，式(13)具有3阶代数精度，且在任意积分区间[a,b]上截断误差为

(15)

为了对比，给出同阶三点和四点Newton-Cotes公式的截断误差分别为

(16)

(17)

由此可见，新四点积分公式具有更小的截断误差系数绝对值，因而数值积分的精度要优于同阶的Newton-Cotes公式。

至此，利用式(13)完全可以采用王永等的思路实现式(2)中对Duhamel积分项的数值积分。

采用4阶显式龙格-库塔法计算预估值yk+i/4(i=1,3,4)如下

(18)

针对文献[6]不足，文中所提出的改进精细积分单步法，在计算量上改进精细积分单步法也只需要进行一次指数矩阵的精细积分即可，其余指数矩阵可通过矩阵乘法获得。

2 算例分析

算例1 二自由度线性动力学方程

(19)

将其整理为式(1)的形式，则有：

(20)

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

不同算法的仿真步长都取Δt=0.2 s，积分时程为10 s，预估公式都采用经典四阶Runge-Kutta法，分别采用本文方法(简记为BM3)、精细积分-微分求积法(简记为DQ4)和精细科茨法(精细积分法加辛普森公式简记为NC2，精细积分法加3阶代数精度的科茨公式简记为NC3)求解结果的相对于解析解的绝对误差的绝对值R(x1)曲线如图1所示。

图1 位移的绝对误差曲线

从图1中可知，在代数精度相同的前提下，本文方法由于具有更小的截断误差系数绝对值，因而数值积分的精度更高。

(21)

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

以椭圆积分得到该问题的解析解为基准，分别采用本文方法(BM3)和精细积分-微分求积法(简记为DQ4)计算幅角v1，计算步长取Δt=0.01 s，相应的计算结果列入表1中。

表1 幅角v1的数值结果对比(Δt=0.01 s)

从表1可知，本文方法的计算精度要明显高于精细积分-微分求积法，且与椭圆积分解完全吻合，这再次验证了本文方法在计算精度上的优势。

算例3 考虑Van der Pol振子方程为

(22)

(23)

表2 Van der Pol方程中x的计算结果对比(Δt=0.1 s)

从表2可知，本文方法的计算结果可以与细分解保持小数点后六位保持一致，计算精度要明显高于同阶的精细积分-微分求积法和预估校正-辛时间子域法[14]。表3是此时本文算法和文献[5]中算法的计算时间对比，从中可知本文算法的计算效率较高，这种优势随着计算步数的增加表现得更为突出。

表3 文献[5]中单步法与本文单步法的计算时间比较

图2 Van der Pol振子位移轨迹(BM3)

图3 Van der Pol振子位移轨迹(TR)

从图2～3可知，对于非线性刚性Van der Pol振子方程，采用本文方法也可以获得比较精确的数值解。这说明虽然本文方法严格上属于显式积分算法，但是相对于更适合求解刚性微分方程的隐式梯形算法依然具有较大的优势。

算例4 考虑平方非线性的二自由度动力方程为

(24)

(25)

取待求变量的初值v1(0)=v2(0)=0.1，v3(0)=v4(0)=0。本文方法的步长分别在Δt=0.1 s和Δt=0.5 s时，细分解和本文方法以及DQ4的数值结果比较如表4所示。

表4 非线性二自由度动力学方程x、y计算结果对比

通过表4可知，本文方法在步长较小时有着与细分解相当的计算精度，DQ4的计算精度与BM3相当，继续将仿真步长扩大五倍，BM3仍能取得与细分解较为接近的结果。

3 结 论

本文是对文献[6]中思路的延续和改进，基于单步块方法和精细积分法得出了如下的结论：

(1) 采用均匀网格的单步块方法构造的数值积分公式即是Newton-Cotes公式；

(2) 采用CGL网格的单步块方法构造的四点数值积分公式比同阶代数精度的Newton-Cotes公式的数值积分精度要高；

(3) 本文提出基于单步块方法的改进精细积分单步法比预估校正-辛时间子域法、文献[5]中的单步法具有更高的计算精度和效率，比文献[6]中的单步法计算精度要高，计算效率理论上两者相当；

(4) 结合文献[6]中方法和改进精细积分单步法可以构成自适应显式精细积分单步法，该方法可以应用于结构地震碰撞反应分析[17-18]。

附录A

Padé逼近是一种有理数逼近，它克服了用多项式逼近大挠度函数效果不理想，而用幂函数(如Taylor级数)逼近函数收敛性太差等缺点。记

(A.1)

Hm为m次多项式构成的集合，则R(m,n)={R∶R=P(x)/Q(x);P(x)∈Hm,Q(x)∈Hn}，分子为m次多项式，分母为n次多项式(除去恒为零的元素)的有理分式的集合。Padé逼近法是从幂级数出发获得有理逼近的一种十分简洁而且非常有效的方法，其基本思想是对于一个给定的形式幂级数，构造一个有理函数，使该有理函数的Taylor展开有尽可能多的项和原来的幂级数相吻合。其具体定义如下

定义：如果有理函数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。

依据Padé逼近的定义，f(x)的(m,n)-Padé逼近也可以理解为由方程

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

(A.2)

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