ABSTRACT

A Highly Accurate Fox-Goodwin time Integration Method (HAFIM) is presented in this paper. HAFIM creates the highly accurate Jacobi matrix by multi-sub-step notion and uses the highly accurate Jacobi matrix for recursive procedure. In the process of constructing the matrix, the exponential matrix operation techniques are employed to reduce the computational efforts and the computer rounding error. The analysis on properties indicates that the proposed method possesses higher amplitude and phase accuracy compared with the Precise time Integration Methods (PIMs). Moreover, the proposed method is very practical given its simplicity.

Keywords:Fox-Goodwin Method, High Accuracy, Time Integration Method, Precise Integration Method

高精度Fox-Goodwin时间积分法

邢修源¹，张慧敏²

¹中国矿业大学(北京)附属中学，北京

²北京航空航天大学，北京

收稿日期：2017年11月22日；录用日期：2017年12月6日；发布日期：2017年12月12日

摘 要

本文建立了高精度Fox-Goodwin时间积分法。这种方法采用等分步的思想，连乘每分步的Jacobi矩阵建立高精度Jacobi矩阵，然后利用高精度Jacobi矩阵进行递推求解。在构造高精度Jacobi矩阵的过程中，用指数矩阵运算技巧来达到降低计算量，用存贮增量矩阵策略来减小舍入误差的影响。性能分析表明，这种方法与精细积分法相比，具备较高的幅值和相位精度，且格式简单，便于应用。

关键词 :Fox-Goodwin方法，高精度，时间积分法，精细积分法

Copyright © 2017 by authors and Hans Publishers Inc.

This work is licensed under the Creative Commons Attribution International License (CC BY).

http://creativecommons.org/licenses/by/4.0/

1. 引言

时间积分法是求解结构动力学响应强有力的方法之一，它的基本思想是将时间域离散，通过迭代格式求出各离散时刻的物理量近似值。经典时间积分法有Houbolt方法 [1] ，Newmark方法 [2] ，Wilson-θ方法 [3] ，Hilber-Hughes-Taylor-α方法 [4] 和Wood-Bossak-Zienkiewicz-α方法 [5] 等，上世纪九十年代以来也出现了诸多新颖的时间方法，例如Zhai方法 [6] ，广义-α方法 [7] ，Bathe方法 [8] 和基于加权残量思想的方法 [9] 等。但这些方法大都只有二阶精度，需要使用非常小的步长才能实现高精度要求，导致大量的计算量。

对于线性定常系统，钟 [10] 提出了精细积分方法(Precise Integration Method, PIM)，以其高精度，高效率等优点得到了广泛应用。PIM的基本思想是对微分方程的解析解进行精确数值计算，对于齐次方程，PIM可以给出计算机精度的解。而对于非齐次方程，需要求解非齐次项产生的Duhamel积分，钟 [11] 和林 [12] 就非齐次项分别为多项式、弦函数、指数函数以及它们的组合形式时，推导了相应的精细积分递推格式，后来又发展了很多种非齐次项的积分方法，例如直接对非齐次项Duhamel积分进行数值积分求解 [13] ，或者增维精细积分方法 [14] 等。当用于求解结构动力学响应时，PIM是一种条件稳定算法 [15] ，它的精度会随着展开项的增多而提高。

本文借助PIM中使用的指数矩阵乘法和存储增量的技术，基于Fox-Goodwin方法 [16] 构造高精度时间积分法Highly Accurate Fox-Goodwin time Integration Method (HAFIM)。Fox-Goodwin方法是一种辛几何算法，对于无阻尼系统可实现四阶精度，并且具有足够的稳定区间。对于收敛的算法来说，当步长趋近于0时，数值解趋近于解析解。HAFIM的思想是将每个时间步长N等分，在每一个分步内使用Fox-Goodwin方法，通过连乘每分步的Jacobi矩阵计算出单步的高精度Jacobi矩阵，利用该Jacobi矩阵进行迭代求解。在构造高精度Jacobi矩阵的过程中，2^m算法和增量算法分别被用来减少计算量和降低舍入误差。与PIM相比，HAFIM具有更简单的矩阵形式和更优秀的数值性能。

2. 算法格式

2.1. 动力学方程

结构动力学方程的一般形式为

$M\ddot{x}+C\dot{x}+Kx=R\left(t\right)$ (1)

式中，M，C，K分别为质量，阻尼和刚度矩阵， $\ddot{x}$ ， $\dot{x}$ ，x分别为加速度，速度和位移向量，R(t)为外部激励。若R(t)可表示为一系列基本函数的组合，且这些函数可以用其本身和一阶和二阶导数线性表出，即

$R\left(t\right)=Tf\left(t\right),f\in \left\{f\left(t\right)|\ddot{f}+B\dot{f}+Df=0,B,D为常矩阵\right\}$ (2)

则可使用增维的方式将外载荷函数转化为待求变量，从而将非齐次方程(1)转化为齐次方程

$\left[\begin{array}{cc}M& O\\ O& I\end{array}\right]\left[\begin{array}{c}\ddot{x}\\ \ddot{f}\end{array}\right]+\left[\begin{array}{cc}C& O\\ O& B\end{array}\right]\left[\begin{array}{c}\dot{x}\\ \dot{f}\end{array}\right]+\left[\begin{array}{cc}K& -T\\ O& D\end{array}\right]\left[\begin{array}{c}x\\ f\end{array}\right]=\left[\begin{array}{c}0\\ 0\end{array}\right]$ (3)

式中，I为单位矩阵，O为零矩阵。下面给出最一般的转化方式。

在有限时间内，任意连续函数都可展开成傅里叶级数的形式，则R(t)可近似表示为

$R\left(t\right)\cong a_0+{\displaystyle \sum _{n=1}^q\left(a_n\cos \left(\frac{2n\text{π}t}{T_0}\right)+b_n\cos \left(\frac{2n\text{π}t}{T_0}\right)\right)}$ (4)

式中，T₀为R(t)的基本周期，a₀，a_n和b_n为傅里叶系数，q为展开项数。若令

$\begin{array}{l}T=\left[\begin{array}{ccccccc}{a}_0& a_1& \cdots & a_q& b_1& \cdots & b_q\end{array}\right]\\ f={\left[\begin{array}{ccccccc}1& \sin \left(\frac{\text{π}t}{T_0}\right)& \cdots & \sin \left(\frac{q\text{π}t}{T_0}\right)& \cos \left(\frac{\text{π}t}{T_0}\right)& \cdots & \cos \left(\frac{q\text{π}t}{T_0}\right)\end{array}\right]}^{\text{T}}\\ B=O\\ D=\left[\begin{array}{ccccccc}0& 0& \cdots & 0& 0& \cdots & 0\\ 0& {\left(\frac{\text{π}t}{T_0}\right)}^2& \cdots & 0& 0& \cdots & 0\\ ⋮& ⋮& \ddots & ⋮& ⋮& \ddots & ⋮\\ 0& 0& \cdots & {\left(\frac{q\text{π}t}{T_0}\right)}^2& 0& \cdots & 0\\ 0& 0& \cdots & 0& {\left(\frac{\text{π}t}{T_0}\right)}^2& \cdots & 0\\ ⋮& ⋮& \ddots & ⋮& ⋮& \ddots & ⋮\\ 0& 0& \cdots & 0& 0& \cdots & {\left(\frac{q\text{π}t}{T_0}\right)}^2\end{array}\right]\end{array}$ (5)

则任意函数均可近似转化为方程(2)的形式。特别地，如果外激励函数本身是二阶齐次常微分方程的解，如常激励，简谐激励，或它们的组合形式等，方程(2)可严格满足。

基于以上讨论，本文仅关注齐次方程的求解，即

$\tilde{M}\ddot{\tilde{x}}+\tilde{C}\dot{\tilde{x}}+\tilde{K}\tilde{x}=0$ (6)

式中

$\tilde{x}=\left[\begin{array}{c}x\\ f\end{array}\right],\tilde{M}=\left[\begin{array}{cc}M& O\\ O& I\end{array}\right],\tilde{C}=\left[\begin{array}{cc}C& O\\ O& B\end{array}\right],\tilde{K}=\left[\begin{array}{cc}K& -T\\ O& D\end{array}\right]$ (7)

2.2. 求解格式

时间积分法求解方程(6)的简明递推格式为

$z_{k+1}=A\left(h\right)z_k$ (8)

式中，A为Jacobi矩阵或称为传递矩阵，z_k为在时间点 $t_k=t_0+kh$ 的状态向量，h为时间步长。本文采用Fox-Goodwin方法，其差分格式为

$\begin{array}{l}\tilde{M}{\ddot{\tilde{x}}}_{k+1}+\tilde{C}{\dot{\tilde{x}}}_{k+1}+\tilde{K}{\tilde{x}}_{k+1}=0\\ {\tilde{x}}_{k+1}={\tilde{x}}_k+h{\dot{\tilde{x}}}_k+h^2\left(\frac{5}{12}{\ddot{\tilde{x}}}_k+\frac{1}{12}{\ddot{\tilde{x}}}_{k+1}\right)\\ {\dot{\tilde{x}}}_{k+1}={\dot{\tilde{x}}}_k+h\left(\frac{1}{2}{\ddot{\tilde{x}}}_k+\frac{1}{2}{\ddot{\tilde{x}}}_{k+1}\right)\end{array}$ (9)

则对应的Jacobi矩阵为

$\begin{array}{l}A\left(h\right)=\left[\begin{array}{cc}I-\frac{h^2}{2}{A}_1+\frac{h^3}{12}{A}_3\left(A_2{A}_1+\frac{h}{2}{A}_1^2\right)& hI-\frac{h^2}{2}{A}_2-\frac{h^3}{12}{A}_3\left(A_1-A_2^2-\frac{h}{2}{A}_1{A}_2\right)\\ -h{A}_1+\frac{h^2}{2}{A}_3\left(A_2{A}_1+\frac{h}{2}{A}_1^2\right)& I-h{A}_2-\frac{h^2}{2}{A}_3\left(A_1-A_2^2-\frac{h}{2}{A}_1{A}_2\right)\end{array}\right]\\ z_k={\left[\begin{array}{cc}{\tilde{x}}_k^{\text{T}}& {\dot{\tilde{x}}}_k^{\text{T}}\end{array}\right]}^{\text{T}},A_1={\tilde{M}}^{-1}\tilde{K},A_2={\tilde{M}}^{-1}\tilde{C},A_3={\left(\tilde{M}+\frac{h}{2}\tilde{C}+\frac{h^2}{12}\tilde{K}\right)}^{-1}\tilde{M}\end{array}$ (10)

HAFIM用多分步的思想，将步长h等分为N份，在每个分步内使用时间积分法，则其递推格式可写为

$z_{k+1}=A_{\text{H}}\left(h\right)z_k$ (11)

式中，A_H为高精度Jacobi矩阵，可由各分步的Jacobi矩阵连乘得到，即

$A_H(h)=\underset{N}{\underbrace{A\left(h_N\right)\cdots A\left(h_N\right)}}=A^N\left(h_N\right)$ (12)

式中，分步长 $h_N=h/N$ 。在构造A_H的过程中，2^m算法和增量算法可分别用来减小计算量和计算机舍入误差。

1) 2^m算法

取 $N={\text{2}}^m$ ，如 $m=\text{2}0$ ， $N=1048576$ ，则仅需m次迭代就可得到A_H，相当于执行下列语句

$\text{for}\left(i=0;i<m;i++\right)A\left(h_N\right)=A\left(h_N\right)\times A\left(h_N\right)$ (13)

2) 增量算法

当h_N很小时，A(h_N)与单位矩阵十分接近，可写为

$A\left(h_N\right)=I+S\left(h_N\right)$ (14)

式中，S(h_N)为增量矩阵

$\begin{array}{l}S\left(h_N\right)=h_N\left[\begin{array}{cc}-\frac{h_N}{2}{A}_1+\frac{h_N^2}{12}{A}_3\left(A_2{A}_1+\frac{h_N}{2}{A}_1^2\right)& I-\frac{h_N}{2}{A}_2-\frac{h_N^2}{12}{A}_3\left(A_1-A_2^2-\frac{h_N}{2}{A}_1{A}_2\right)\\ -A_1+\frac{h_N}{2}{A}_3\left(A_2{A}_1+\frac{h_N}{2}{A}_1^2\right)& -A_2-\frac{h_N}{2}{A}_3\left(A_1-A_2^2-\frac{h_N}{2}{A}_1{A}_2\right)\end{array}\right]\\ A_1={\tilde{M}}^{-1}\tilde{K},A_2={\tilde{M}}^{-1}\tilde{C},A_3={\left(\tilde{M}+\frac{h_N}{2}\tilde{C}+\frac{h_N^2}{12}\tilde{K}\right)}^{-1}\tilde{M}\end{array}$ (15)

可以看出，当h_N趋近于0时，S(h_N)接近于零矩阵。由于一般台式计算机仅有16位有效数字，S(h_N)的精度可能会在迭代过程中丧失殆尽。为避免精度损失，在构造A_H的过程中，要单独存储S(h_N)，而不是A(h_N)。又因为

$\left(I+S\left(h_N\right)\right)\times \left(I+S\left(h_N\right)\right)=I+2S\left(h_N\right)+S\left(h_N\right)\times S\left(h_N\right)$ (16)

则需执行的语句变为

$\text{for}\left(i=0;i<m;i++\right)S\left(h_N\right)=2S\left(h_N\right)+S\left(h_N\right)\times S\left(h_N\right)$ (17)

经过m次乘法后，S不再是一个很小的矩阵，此时A_H可由下式得到

$A_{\text{H}}\left(h\right)=I+S\left(h_N\right)$ (18)

综上所述，表1给出了HAFIM求解线性定常系统的执行步骤。在计算过程中，HAFIM与PIM有两处不同：首先，HAFIM直接给出了增量矩阵S(h_N)的形式，而PIM的增量矩阵需由系数矩阵进一步计算得到；其次，HAFIM无需对结果进行变换，可直接求得位移、速度和加速度，而PIM的速度和加速度需由进一步变换得到，所以说HAFIM计算更为简便，所需时间成本也更低。

实际上，任意一种时间积分法都可按照同样的方式来构造相应的高精度方法，本文选取Fox-Goodwin方法是由于其具备良好的精度和稳定性，下面我们将说明这种方法的性能优势。

表1. 高精度Fox-Goodwin时间积分方法的执行步骤

3. 算法性能

衡量算法性能优劣的指标主要有稳定性和精度两个方面。对于线性定常系统，由于模态正交性，可以通过对单自由度系统问题的分析来说明算法的特性，即振动方程为

$\ddot{x}+2\xi \omega +{\omega }^{2}x=0$ (19)

式中，ξ和ω分别为阻尼率和固有频率。为避免参数N的影响，这里仅研究Fox-Goodwin方法本身的性质，并与梯形法则(Trapezoidal-Rule, TR)，中心差分法(Central-Difference-Method, CDM)和PIM进行对比。TR和CDM是结构动力学中常用的时间积分法，对无阻尼系统不具有数值阻尼，精度较高。

3.1. 稳定性

稳定性要求算法的谱半径 $\rho \left(A\right)=\text{max}\left\{{\lambda }_{\text{1}},{\lambda }_{\text{2}}\right\}\le \text{1}$ 。对于单自由度方程(19)，Fox-Goodwin方法的Jacobi矩阵为

$A_{\text{F}}=\frac{1}{{\tau }^2+12\xi \tau +12}\left[\begin{array}{cc}-4\xi {\tau }^3-5{\tau }^2+12\xi \tau +12& -8{\xi }^2{\tau }^2+12\\ 2{\tau }^4-12{\tau }^2& 4\xi {\tau }^3-5{\tau }^2-12\xi \tau +12\end{array}\right],z_k={\left[\begin{array}{cc}{x}_k& h{\dot{x}}_k\end{array}\right]}^{\text{T}}$ (20)

式中， $\tau =\omega h$ 。类似地，TR和CDM的Jacobi矩阵分别如下

$A_{\text{T}}=\frac{1}{{\tau }^2+4\xi \tau +4}\left[\begin{array}{cc}-{\tau }^2+4\xi \tau +4& 4\\ -4{\tau }^2& -{\tau }^2-4\xi \tau +4\end{array}\right],z_k={\left[\begin{array}{cc}{x}_k& h{\dot{x}}_k\end{array}\right]}^{\text{T}}$ (21)

$A_{\text{C}}=\frac{1}{4\xi \tau +4}\left[\begin{array}{cc}-2\xi {\tau }^3-2{\tau }^2+4\xi \tau +4& -4{\xi }^2{\tau }^2+4\\ {\tau }^4-4{\tau }^2& 2\xi {\tau }^3-2{\tau }^2-4\xi \tau +4\end{array}\right],z_k={\left[\begin{array}{cc}{x}_k& h{\dot{x}}_k\end{array}\right]}^{\text{T}}$ (22)

而PIM的性质随着展开项数的增多而不同，这里给出最常用的两种情况：

1) L = 3

$A_3=I+hH+\frac{h^2}{2}{H}^2+\frac{h^3}{6}{H}^3$ (23)

2) L = 4

$A_4=I+hH+\frac{h^2}{2}{H}^2+\frac{h^3}{6}{H}^3+\frac{h^4}{24}{H}^4$ (24)

式中L表示截断阶数，而

$H=\left[\begin{array}{cc}-\xi \omega & 1\\ {\xi }^2{\omega }^2-{\omega }^2& -\xi \omega \end{array}\right],z_k={\left[\begin{array}{cc}{x}_k& {\dot{x}}_k+\xi \omega x_k\end{array}\right]}^{\text{T}}$ (25)

图1和图2分别给出了无阻尼情况和阻尼率为0.1时这几种方法的谱半径曲线。可以看出，TR为无条件稳定算法，其余方法均为条件稳定方法，Fox-Goodwin方法的稳定区间介于两种精细积分法之间。为方便应用，表2给出了这三种方法的稳定极限τ_cr，当 $\tau \le {\tau }_{\text{cr}}$ 时，算法可以给出稳定的结果。

3.2. 精度

考虑单自由度系统(19)，该问题的解析解为

$x\left(t\right)={\text{e}}^{-\xi \omega t}\left(c_1\cos {\omega }_{\text{d}}t+c_2\sin {\omega }_{\text{d}}t\right)$ (26)

式中，c₁和c₂为由初始条件确定的常数， ${\omega }_{\text{d}}=\omega \sqrt{1-{\xi }^2}$ 。而数值解的形式为

$x_k={\text{e}}^{-\overline{\xi }\overline{\omega }{t}_k}\left({\overline{c}}_1\cos {\overline{\omega }}_{\text{d}}{t}_k+{\overline{c}}_2\sin {\overline{\omega }}_{\text{d}}{t}_k\right)$ (27)

式中

图1. 谱半径曲线( ξ = 0)

图2. 谱半径曲线( $\xi$ = 0.1)

表2. 算法的稳定极限

$\begin{array}{l}{\lambda }_{1,2}=A\pm B\text{i}={\text{e}}^{\left(-\overline{\xi }\overline{\omega }\pm \text{i}{\overline{\omega }}_{\text{d}}\right)h}\\ {\overline{\omega }}_d=\frac{\arctan \left(B/A\right)}{h}\\ \overline{\omega }=\frac{{\overline{\omega }}_{\text{d}}}{\sqrt{1-{\overline{\xi }}^2}}\\ \overline{\xi }=-\frac{\ln \left(A^2+B^2\right)}{2\overline{\omega }h}\end{array}$ (28)

在稳定区间内，λ_1,2为一对共轭复根。将数值解的形式与解析解对比，我们可以用 $\overline{\xi }$ 来表示数值耗散的程度，称为幅值衰减率。而数值弥散一般用周期延长率来表示，定义为

$\eta =\frac{\overline{T}-T}{T},\overline{T}=\frac{2\text{π}}{{\overline{\omega }}_{\text{d}}},T=\frac{2\text{π}}{{\omega }_{\text{d}}}$ (29)

图3和图4分别给出了这几种格式的幅值衰减率和周期延长率曲线，结果表明无阻尼情况下，PIM的幅值精度较差，TR和CDM方法的相位精度较差，而Fox-Goodwin方法的幅值和相位精度皆优于其它方法，因此本文选用Fox-Goodwin方法来构造高精度时间积分法。特别地，由于Fox-Goodwin方法的辛几何特性，其幅值衰减率为零，这意味着幅值误差不会随着步数的增多而累积，对长期仿真十分有利。

4. 数值仿真

考虑三自由度无阻尼质点-弹簧系统，如图5所示。弹簧的刚度系统均为k = 1 N/m²，质点的质量均为m = 1 kg，质点1处受到简谐激励R = 10sin5t N作用，系统的运动方程为

$\left[\begin{array}{ccc}1& 0& 0\\ 0& 1& 0\\ 0& 0& 1\end{array}\right]\left[\begin{array}{c}{\ddot{x}}_1\\ {\ddot{x}}_2\\ {\ddot{x}}_3\end{array}\right]+\left[\begin{array}{ccc}2& -1& 0\\ -1& 2& -1\\ 0& -1& 2\end{array}\right]\left[\begin{array}{c}{x}_1\\ x_2\\ x_3\end{array}\right]=\left[\begin{array}{c}10\sin 5t\\ 0\\ 0\end{array}\right]$ (30)

令 $x_{\text{4}}=\text{sin5}t$ ，则方程(28)可转化为齐次方程

$\left[\begin{array}{cccc}1& 0& 0& 0\\ 0& 1& 0& 0\\ 0& 0& 1& 0\\ 0& 0& 0& 1\end{array}\right]\left[\begin{array}{c}{\ddot{x}}_1\\ {\ddot{x}}_2\\ {\ddot{x}}_3\\ {\ddot{x}}_4\end{array}\right]+\left[\begin{array}{cccc}2& -1& 0& -10\\ -1& 2& -1& 0\\ 0& -1& 2& 0\\ 0& 0& 0& 25\end{array}\right]\left[\begin{array}{c}{x}_1\\ x_2\\ x_3\\ x_4\end{array}\right]=\left[\begin{array}{c}0\\ 0\\ 0\\ 0\end{array}\right]$ (31)

初始条件为

$x_0={\left[\begin{array}{cccc}0& 0& 0& 0\end{array}\right]}^{\text{T}},{\dot{x}}_0={\left[\begin{array}{cccc}0& 0& 0& 5\end{array}\right]}^{\text{T}}$ (32)

图3. 幅值衰减率曲线( ξ = 0)

图4. 周期延长率曲线( ξ = 0)

在本例中，由于激励为简谐函数，齐次方程(29)与原方程(28)的解析解完全相同。

简谐激励周期为 $T=\text{2π}/\omega =\text{1}.\text{2s}$ ，因此采用的时间步长选为 $h=T/10=0.\text{1s}$ ，令 $m=\text{5}$ ，应用HAFIM和PIM (L = 3, L = 4)计算了该系统在[0, 40 s]内的位移，速度和加速度响应。质点1的结果如图6~8所示，其中精确解由模态叠加法得到。可以看出，由于HAFIM的保辛特性，其幅值误差不累积；而经过一段

图5. 质点–弹簧系统

图6. 质点1的位移响应和数值方法的绝对误差

图7. 质点1的速度响应和数值方法的绝对误差

图8. 质点1的加速度响应和数值方法的绝对误差

时间后，PIM的幅值耗散会大大降低结果的精度，尤其是位移和加速度响应更加明显，因此HAFIM对长期仿真十分有利。此外，虽然PIM的位移精度较高，但其速度精度较差，即使是截断到第4阶(L = 4)，其速度误差也远大于HAFIM给出的结果。

5. 结论

本文构造了高精度Fox-Goodwin时间积分法。它的基本思想是将每个时间步等分为N份，连乘每个分步的Jacobi矩阵构造高精度Jacobi矩阵，然后使用高精度Jacobi矩阵进行递推求解。在构造高精度Jacobi矩阵的过程中，2^m算法和增量算法被用来提升计算效率，减小舍入误差。

相比较于精细积分法来说，本文提出的方法结合了时间积分法本身的优势，使得其在计算量和性能上得到优化。首先，这种方法的增量矩阵无需进一步计算，可采用文中给出的形式，并可直接得到位移和速度的数值解；其次，Fox-Goodwin方法的辛几何特性使得该方法具备较高的幅值精度，且应用于无阻尼系统时是一种四阶格式；最后，数值结果证实了高精度Fox-Goodwin时间积分法的精度优势。

致谢

感谢国家自然科学基金资助项目(批准号为11372021和11672019)的资助。

文章引用

邢修源,张慧敏. 高精度Fox-Goodwin时间积分法
Highly Accurate Fox-Goodwin Time Integration Method[J]. 力学研究, 2017, 06(04): 161-172. http://dx.doi.org/10.12677/IJM.2017.64017

参考文献 (References)