Page 326 - 《振动工程学报》2025年第11期

2784 振 动 工 程 学 报 第 38 卷

行为。为了降低其超调行为，计算 c 2 = 0可得：  b
A = 1
 3
√ 


( ρ +8+ρ ∞ )(1−ρ ∞ )−4  (2α 1 + p+2pα 6 )Ω +2−6p
2
2


α 6 = ∞ （17）  b b
A =

4(ρ ∞ +1)(ρ ∞ −2)  2 2p





本文提出了一类自启动单解隐式算法，它实现 −(2α 1 −1+2pα 6 )Ω −2+3p （18）
2

 b b
 A =

了一致二阶精度、可控数值高频耗散、无条件稳定、  1 p





 2
可接受的超调行为、最优的振幅与相位误差主项。  b (2α 1 − p+2pα 6 )Ω +2−2p

b

 A =

 0
上述过程也可以较容易地扩展到隐式算法的三参数 2p
显 式 算 法 实 现 分 岔 点 处 可 控 数 值 耗 散 的 条
版本 ( λ , j = 1, 2, 3)，故此处不再赘述。
∞
j

件 [2,12] 为：
b b 2 b 2 （19）
3 显 式 算 法 的 优 化 A = 2ρ b +ρ s , A = 2ρ b ρ s +ρ , A = ρ ρ s
1
0
b
2
b
式中， 0 ⩽ |ρ s | ⩽ ρ b ⩽ 1为两个用户指定参数，用于控制
显式算法在稳定域内的数值耗散量。将式 (18) 代入
当 α 2 = 0时，式 (6) 退化为一致二阶精度的显式
式 (19) 并求解可得：
算 法 。 此 时， 显 式 算 法 的 数 值 放 大 矩 阵 的 不 变 量
√
A i (i = 0, 1, 2, 3)在分岔点 ( Ω = Ω b ) 处化简为： Ω b = (1+ρ b )(2+ρ s −ρ b ρ s ) （20）

2+ρ s −ρ b ρ s

 p =


 (1+ρ b )(1+ρ s )

 （21）
 2 2
 2(1+ρ b ) (1+ρ s )(ρ b ρ s −ρ s −2)α 1 +(1−ρ b )ρ s (ρ ρ s +ρ s +4)+2(ρ b +3)
 b

α 6 =

 2
2(1+ρ b )(2+ρ s −ρ b ρ s )
目前，显式算法的 α 1 和 α 4 仍是未知的。文献 [20-21] 位误差由式 (12) 表出，其中误差主项分别为：
[ ]
通过泰勒级数截断选取 α 1 = p /2和 α 4 = p/2，但该组 12pα 4 −6p+1 2 3
2
δ 2 = p(α 4 +α 6 −1)+α 1 − ξ ξω （22）
取值并未实现误差的最小化。显式算法的振幅与相 3
4
2
4(18pα 4 +6pα 6 −12p+6α 1 +1)ξ −12α 1 −8(12pα 4 −6p+1)ξ −12pα 6 +6p+1
ε 2 = √ ω 3 （23）
24 1−ξ 2
(4)
2
3
(3)
2
式中， p和 α 6 由式 (21) 给出。由于参数 α 4 决定了二阶 τ 1 = (ζ 1 u +ζ 2 ξωu +ζ 3 ω ¨u n )∆t +O(∆t ) （25）
n
n
(k)
微分方程中预测速度格式的显隐式属性，故需要分 式中， ¨ u n 为单自由度系统的加速度； u 表示位移的
n
以下两类情形： α 4 = 0和 α 4 , 0。 k阶导数在 时刻的取值且系数 ζ j (j = 1, 2, 3)为：
t n

 1+6p(1−2α 6 )−6p 2

 ζ 1 =


3.1 完全显式算法  12p




 2

 1−6pα 6 −3p
 （26）
当 α 4 = 0时，显式算法采用了速度的显式处理， ζ 2 = 3p






退化为完全显式算法：  2α 1 − p 2


 ζ 3 =


p α 1 0 p 0 2p
1 1 1 1 1 速度和加速度的局部截断误差亦有类似的表达
−α 6 α 6 1− 1−
2 2p 2p p p 形式，但其系数值独立于参数 α 1 ，因此只能考虑位移
（24）
的局部截断误差确定参数 α 1 。
式中， p和 α 6 由式 (21) 给出。
(1) 当通过计算 ζ 3 = 0获得参数 α 1 时，其值退化为
该完全显式算法的振幅与相位误差主项 ( δ 2 和
文献 [20] 中的设计思想：
ε 2 ) 独立于参数 α 1 。换句话说，本文无法通过优化振 p 2 (2+ρ s −ρ b ρ s ) 2
α 1 = = （27）
幅或相位误差主项确定参数 α 1 。 2 2(1+ρ b ) (1+ρ s ) 2
2
计算显式算法的位移局部截断误差 [12] 为： (2) 当通过计算 ζ 2 = 0获得参数 α 1 时，其值为：
3
3
3
2
2
2
3
2
2
(5ρ +ρ s −1)ρ s ρ −5ρ −29ρ −(13ρ +21ρ −3ρ s −2)ρ −61ρ s +(13ρ +57ρ +51ρ s −2)ρ b −40
s
s
s
s
s
s
s
α 1 = b b （28）
2
2
6(1+ρ s ) (1+ρ b ) (ρ b ρ s −ρ s −2)
(3) 当通过计算 ζ 1 = 0获得参数 α 1 时，其值为：
3
3
3
2
2
3
2
2
2
(17ρ +10ρ s −1)ρ s ρ −5ρ −32ρ −(31ρ +72ρ +21ρ s −2)ρ −79ρ s +(19ρ +102ρ +117ρ s +16)ρ b −58
s
s
s
s
s
s
s
α 1 = b b （29）
2
2
12(1+ρ s ) (1+ρ b ) (ρ b ρ s −ρ s −2)
(4) 当定义误差总量 Γ = ζ +ζ +ζ 且通过求解 ∂Γ/∂α 1 = 0最小化 Γ时，参数 α 1 计算为：
2
2
2
3
2
1
3
3
2
2
3
2
2
2
3(7ρ +2ρ s −1)ρ s ρ −17ρ −100ρ −(51ρ +92ρ −ρ s −6)ρ −213ρ s +(47ρ +210ρ s +119)ρ s ρ b −142
s
s
s
s
s
s
α 1 = b b （30）
2
2
24(1+ρ s ) (1+ρ b ) (ρ b ρ s −ρ s −2)