第 62 卷第 5 期2023 年9 月
Vol.62 No.5
Sept.2023中山大学学报（自然科学版）（中英文）
ACTA SCIENTIARUM NATURALIUM UNIVERSITATIS SUNYATSENI
基于Tikhonov正则化改进的IHB法求解
Mathieu-Duffing系统多重解*
王德亮1，2，刘济科1，刘广1，2
1. 中山大学航空航天学院，广东深圳 518107
2. 深圳市智能微小卫星星座技术与应用重点实验室，广东深圳 518107
摘要：增量谐波平衡法（IHB法）是研究强非线性振动系统的一种半数值半解析方法，然而已有研究表明，在求解含多重解的系统时该方法的收敛性强烈地依赖于初值的选择。Tikhonov正则化常被用于优化问题中来解决可能出现的病态问题。文章通过在原始的IHB法中引入Tikhonov正则化，提出一种改进的IHB法
（TIHB法）来求解具有多重解的Mathieu-Duffing系统。结果表明，改进的TIHB法可以快速、高效地获得系统的多个稳定或不稳定解，且算法的收敛性能要远远优于原始的IHB法。
关键词：非线性振动；IHB法；Tikhonov正则化；多重解
中图分类号：V21 文献标志码：A 文章编号：2097 - 0137（2023）05 - 0078 - 07
Multiple solutions of the Mathieu-Duffing system obtained
by the improved IHB method based on Tikhonov regularization
WANG Deliang1，2， LIU Jike1， LIU Guang1，2
1. School of Aeronautics and Astronautics，Sun Yat-sen University， Shenzhen 518107， China
2. Shenzhen Key Laboratory of Intelligent Microsatellite Constellation， Shenzhen 518107， China
Abstract：The incremental harmonic balance method （IHB method） is a semi-numerical and semi-ana‐lytical method for strongly nonlinear dynamic systems. However， previous studies have shown that the convergence performance of the original IHB method in solving systems with multiple solution
s strong‐ly depends on the selection of initial values. The Tikhonov regularization is often used in optimization problems to solve potential ill-posed problems. In this paper， by incorporating the Tikhonov regulariza‐tion into the original IHB method， an improved IHB method （TIHB method） is proposed to obtain the multiple solutions of the Mathieu-Duffing system. The results show that the improved TIHB method can obtain the stable and unstable solutions of the Mathieu-Duffing system quickly and efficiently， and the convergence performance of the TIHB method is much better than the original IHB method.
Key words：nonlinear vibration； IHB method； Tikhonov regularization； multiple solution
现实中的各种振动系统都含有非线性因素（陈予恕，1992；陈树辉，2007；Amabili，2008）。然而，早期的数学工具和振动理论通常会忽略掉系统的非线性因素而当作线性系统来求解（Devore，1998）。
DOI：10.13471/jki.acta.snus.2023A001
*收稿日期：2023 − 01 − 13 录用日期：2023 − 02 − 27 网络首发日期：2023 − 07 − 13基金项目：国家自然科学基金（11972380）；广东省基础与应用基础研究基金（2021A1515110750，2023A1515010028）；深圳市科技计划（ZDSYS20210623091808026）
作者简介：王德亮（1998年生），男；研究方向：非线性振动；E-mail：*********************.edu
通信作者：刘广（1992年生），男；研究方向：非线性振动、参数识别等；E-mail：****************.edu
第 5 期王德亮，等：基于Tikhonov正则化改进的IHB法求解Mathieu-Duffing系统多重解
某些时候线性近似可能会超过容许误差，甚至导致定性上的错误（Dowell et al.，1988）。随着求解技术尤其是计算技术的发展，非线性振动理论开始逐步完善，其研究方法也日益完善。非线性振动研究可以分为定性和定量两种方法，本文主要关注定量方法。定量方法主要有数值解法和解析方法，此外还有二者相结合的半数值半解析方法。数值方法，如Runge-Kutta法（徐加初等，2008）、Newmark-β法（刘广等，2017）和精细积分法（富明慧等，2008）等。半解析方法有时域最小残值法TMRM（Liu et al.，2021a；Liu et al.，2021b；Liu et al.，2022；Zheng et al.，2022a），同伦分析HAM（Liao et al.，2004）和增量谐波平衡法IHB（张丹伟等，2018；赵倩等，2018）。
IHB法最早由Lau和Cheung et al.于20世纪80年代初提出（Ferri，1986），该方法在应用时不受小参数的限制，是一种适用于弱非线性系统和强非线性系统的半数值半解析方法（Cheung et al.，1990）。IHB法具有概念直观、易于应用的特点，广泛地应用在非线性振动的各种研究中。近30年来，涌现了大量的基于IHB法或改进IHB法的研究。即便如此，原始的IHB法仍有以下一些局限：一方面，IHB法在本质上是谐波平衡法和Newton-Raphson迭代法的结合，其收敛性强烈地依赖于迭代初值的选取（黄建亮
等，2021）；另一方面，对于参数激励的非线性振动系统，可能同时具有多个周期解或极限环解，其中还包含稳定解和不稳定解。在IHB法中，由于线性增量方程的雅可比矩阵可能是病态的，如果参数的初始猜测值远离精确解，计算很可能不收敛（Liu et al.，2018；钱征文等，2011）。而系统同时存在多重解时，IHB法可能难以获得全部的半解析解（Zheng et al.， 2022b）。
在后续的基于IHB法的振动分析中，近些年发展了一些选取迭代初值的思路，比如：结合数值方法获得数值响应，通过快速傅里叶变换（FFT）获得系统的近似频率和振幅，从而确定迭代初值（Ju et al.， 2020）。但这种思路对不稳定解无能为力。另外，可通过求解系统的近似线性系统的级数解作为迭代初值（Chen et al.， 2009）。但某些系统具有丰富的非线性特性，低阶近似解无法准确
揭示这些特性，从而造成方法失效。
本文拟通过在IHB法的线性增量方程中引入Tikhonov正则化，来提高IHB法的收敛性能。并通过改进的IHB法，求解了具有多重解的非线性参数激励Mathieu-Duffing系统。数值结果表明，改进的TIHB不仅可以快速获得高精度的半解析解，而且能通过选取不同的迭代初值获得Mathieu-Duffing系统不同的稳定或不稳定周期解。
1 基于Tikhonov正则化的TIHB法
1.1　IHB法
考虑如下形式的非线性Mathieu-Duffing系统（Shen et al.，2008）：
x+2ζx-(1+αsin2ωt)x+βx3=0 .（1）通过引入时间变换τ=ωt，方程（1）变为
ω2x″+2ζωx′-(1+αsin2τ)x+βx3=0，（2）其中上标“′”表示对无量纲时间τ的导数。对于系统的稳态周期响应，其半解析解为
x=a0+∑k=1∞[]
a k cos kτ+
b k sin kτ .（3）
接下来，通过IHB法求解方程（3）所示的半解析解。IHB法的第一步是增量过程。假设x0为方程（1）的精确解，通过引入增量Δx和Δω，可以得到其邻近点x=x0+Δx和ω=ω0+Δω，代入式（2）得
ω20x″0+2ζω0x′-(1+αsin2τ)x0+βx30+ω20Δx″
+2ζω0Δx′-(1+αsin2τ)Δx+βx20Δx
+2(ω0x″0+ζx′0)Δω+ο(Δx，Δω)=0.
（4）引入方程
R(x0)=ω20x″0+2ζω0x′-(1+αsin2τ)x0+βx30 .（5）方程（4）可简化为
R(x0)+ω20Δx″+2ζω0Δx′
-(1+αsin2τ)Δx+βx20Δx
+2(ω0x″0+ζx′0)Δω+ο(Δx，Δω)=0，(6)
其中R(x
)表示原始方程的残差，ο(Δx，Δω)表示关于增量Δx和Δω的高阶小量。仅保留关于Δx和
79
第 62 卷
中山大学学报（自然科学版）（中英文）Δω的一阶小量，可得增量方程
R (x 0)+ω20Δx ″+2ζω0Δx ′-(1+αsin 2τ)Δx
+βx 20Δx +2(ω0x ″0+ζx ′0)Δω=0 .          (7)
截取Fourier 级数的前N 项作为系统的近似解，得
ì
íîïïïïïïx 0≈x N =a 0+∑k =1N
[]a
k
cos kτ+b k sin kτ，Δx (t )=Δa 0+∑k =1
N
[]a
k
cos kτ+b k sin kτ.
（8）
将方程（8）代入方程（6）中，并引入Galerkin 平
均过程
∫
2π[ω20Δx ″+2ζω0Δx ′-(1+αsin 2τ+3βx 20)Δx
]+2(ω0x ″0+ζx ′0)Δωδ(Δx )d x  =-
∫
2πR (x 0)δ(Δ(x ))d τ，（9）
其中
ìíî
ï
ïïï
ïa =[]a 0，a 1，b 1，⋯，a k ，b k ，Δa =[]Δa 0，Δa 1，Δb 1，⋯，Δa k ，Δb k ，C S =[1，cos τ，⋯，cos Nτ，sin τ，⋯，sin Nτ].（10）
积分式（9）并化简为代数方程组，即可得到关于Δa 、Δω的矩阵形式的增量方程
J a Δa =R
ˉ-J ωΔω，（11）
其中J a 和J ω为方程关于谐波系数和频率的雅可比矩阵，R ˉ为残差（5）积分之后的结果。方程（1）所示
的Mathieu-Duffing 系统为已知周期的参数激励系统，其稳态周期解的频率为ω，因此无需求解。则式（9）可化为J a Δa =R
ˉ，（12）
其中
ìíî
ïïï
ï
ïïï
ïïïï
ïJ a =∫
02π
[ω20C T S C ″S +2ζω0C T
S C ′S ]-C T S (1+αsin 2τ+3β(C S a )2
)C S d τ，R ˉ=-a ∫02πω20C T S C ″S +2ζω0C T S C ′S -C T S (1+αsin 2τ)C S d τ-a ∫
02πβC T S (C S a )2
C S d τ .（13）通过选取合适的迭代初值，并通过方程（12）求解谐波系数的迭代更新量Δa ，然后根据以下方程进行更新
a k =a k -1+Δa ，
（14）直到
Δa <c ，
（15）
则终止迭代。其中c 为预先给定的收敛误差。1.2　引入Tikhonov 正则化的改进IHB 法
在某些的参数条件下，方程（11）中的Jacobi 矩阵J a 和J ω可能有较大的条件数，即谐波系数和频率的Jacobi 矩阵可能是病态的。此时，如果谐波系数的初值选择不当，上述原始的IHB 法的迭代可能不收敛。如何确定IHB 法的迭代初值一直是一个难题，常用的方法可以首先获取系统的稳定数值解，然后通过对数值解进行快速傅里叶变换（FFT ）来获得系统的迭代初值。这个过程无疑是低效且繁琐的，此外对于系统的不稳定解，这种方式还会失效。Tikhonov 正则化经常被用于反问题中来解决可能出现的病态问题，可尝试在方程（12）中引入Tikhonov 正则化，则方程（12）变为Δa λ=arg min Δa
{}
h (a )≔  R
ˉ-J a Δa 2+λ  Δa 2=Δa T (J T a J a +λI )Δa -2Δa T J T a R
ˉ+R ˉT R ˉ .（16）进一步得
Δa λ=[]J T a J a +λI -1J T a R
ˉ .（17）
方程（16）中，  ⋅2
表示2-范数，I 为单位矩阵，λ为正则化参数。从方程（17）中可以看出，不同
的正则化参数λ会获得不同的迭代更新量Δa λ，
一个合适的正则化参数可以平衡计算量和精度之间的矛盾。
常用的确定正则化参数λ的方法有L 曲线法（Ju et al.，2020）。L 曲线法的思想是，选取不同的正则化参数来计算获得方程（16）中  Δa λ2
与  R
ˉ-J a
Δa λ
2
，并且分别以  Δa λ
2
和  R
ˉ-J a Δa λ2作为横纵坐标绘制L 曲线，一般认为L 曲线最大曲
率点对应的正则化参数为最优参数λOPT  .令ρ=log  R
ˉ-J a Δa λ2，η=log  Δa λ2
，则曲率κ的表80
第 5 期王德亮，等：基于Tikhonov正则化改进的IHB法求解Mathieu-Duffing系统多重解达式为
κ=2ρ′η″-ρ″η′
(ρ′2+η′2)3/2，（18）
其中“′”表示对λ的导数。上述通过引入Tik‐
honov正则化改进的IHB法被称为TIHB法。
2 数值算例
2.1　TIHB法求解Mathieu-Duffing系统
对于方程（1）所示的非线性参数激励Mathieu-
Duffing系统，其近似周期解可以通过方程（8）所示
的截断Fourier级数表示，即
x(τ)≈a0+∑k=1N[]
a k cos kτ+
b k sin kτ .（19）
当稳态周期解的频率为参激频率的一半，即
响应周期为激励周期的2倍，此时的稳态周期解称
之为2-周期解。Shen et al.（2008）的研究表明，当正则化英文
选取参数ζ=0.125，β=1，ω=2，3.6<α
<5.264 566时，系统还同时存在1-周期解。上述研
究表明，随着参数的变化，该非线性参数激励Ma‐
thieu-Duffing系统同时存在多个周期解，且周期解
中还有稳定解与不稳定解。
现有的方法在求解这些周期解时存在各种困难。对于数值方法，其无法求得不稳定解；而原始的IHB法由于其收敛性强烈依赖于谐波初值的选择，因此在追踪多个不同解时，需要较为精准地给出不同的迭代初值，这无疑是一项困难的工作。接下来将尝试通过本文提出的TIHB法来解决上述多解问题。进一步选
取α=4，N=40.
本文所提的TIHB法和4阶Runge-Kutta（R-K）法获得的2-周期对应的相图，如图1所示。从图中可以看到，TIHB法所得半解析解和R-K法获得的结果符合得非常好，且即便考虑谐波阶数N=40时，仅仅需要40步迭代就可以收敛。图2给出了TIHB法获得的半解析解对应的残差，也即方程（5）对应的系统残差。从图中可以看到，一个周期内的系统残差在1011量级，这表明TIHB法对应的半解析解具有非常高的精度。
2.2　TIHB法的收敛性能
为了进一步验证本文所提的TIHB法相较于原始IHB法的收敛性能，对方程（19）所示的半解析解的一阶谐波系数赋初值，剩余的谐波系数初值全部设定为0。例如，令一阶谐波系数a1和b1在参数区域[-2.5，2.5]内遍历取值，然后分别通过原始IHB法和本文所提的TIHB法来计算结果。最终的收敛情况如图3所示，其中空心圆表示该迭代初值收敛到系统的“0”解，实心点则表示收敛周期解，空白点则表示该初值不收敛。上述收敛的迭代初值都在40步之内使系统的残差降到10-11以内；空白处的点则是经过500次迭代仍旧无法将残差降到10-10以下的。对比图3（a）和图3（b）可以发现，本文提出的TIHB法在此区域内每一个点都是可以收敛的，但是原始IHB法只是在部分点收敛。上述结果表明，Tikhonov正则化的引入，确实能够显著地增大传统IHB法的收敛域。
结合前文的计算以及对应的收敛性分析，可以发现：相比于传统的IHB法，TIHB法在不损失
精度的前提下，其收敛性得到了极大的改善，或图1　α=4时，2-周期响应的相平面图
Fig.1　Phase diagram of the 2-period response with α=
4
图2　α=4时，2-周期响应对应的残差图
Fig.2　The residual diagram corresponding to
2-periodic response with α=4
81
第 62 卷
中山大学学报（自然科学版）（中英文）者说其收敛域得到了显著的增大。事实上，在当前参数下，系统除了2-周期周期解还同时存在2个1-周期周期解。通过引入约束a 2k -1=b 2k -1=0到方
程（12）中，就可以进一步求得系统的1-周期解。图4为系统的1-周期解sol1和sol2，其中周期解sol1为稳定周期解，周期解sol2为不稳定周期解。
同样地，令一阶谐波系数a 1和b 1在[-2.5，
2.5]内遍历取值，得到原始IHB 法和TIHB 法关于1-周期解的收敛结果，如图5所示。此时系统的收敛情况分为3种，分别是收敛到稳定周期解sol1、不稳
定周期解sol2和“0解”。即便此时的情况比前文中2-周期周期解更复杂，仍旧可以看到TIHB 法各解的收敛域得到了显著的扩大。从图5（a ）中还可以发现，原始IHB 法收敛域中，每个迭代初值似乎是随机收敛到某个解的。
保持Mathieu-Duffing 系统中其他参数不变，取α=5.1。在当前参数设定下，此时系统存在5个周期解，其中两个是不稳定周期解，3个是稳定周期解。通过TIHB 法并结合初值遍历过程，可以依次获得这5个周期解。TIHB 法获得的半解析解对应的相图，如图6所示。其中，周期解sol2和sol5为不稳定解，用虚线表示；其余解为稳定周期解，用实线表示。
进一步，通过原始IHB 法和TIHB 法对α=5.1
的系统进行求解，结果如表1所示。表中的周期解sol1-sol5和图6中的结果相对应。表1显示，原始IHB 法的收敛效果远不如TIHB 法。原始IHB 法在所给的步长和区域内，无法迭代获得解sol1和解sol5，只得到了3个稳定周期解和一个“0”解；且各解的收敛域非常小，在2 601
个初值点域上，仅
图3　原始IHB 法和TIHB 法的收敛性Fig.3　The convergence of the original IHB method
and the TIHB method
图5　原始IHB 法和TIHB 法的收敛性Fig.5　The convergence of the original IHB method
and the TIHB method
图4　系统的1-周期解对应的相图Fig.4　Phase diagram of the 1-period response
82

版权声明：本站内容均来自互联网，仅供演示用，请勿用于商业和其他非法用途。如果侵犯了您的权益请与我们联系QQ:729038198，我们将在24小时内删除。