RECONSTRUCTION HARMONIC BALANCE METHOD AND ITS APPLICATION IN SOLVING COMPLEX NONLINEAR DYNAMICAL SYSTEMS
-
摘要: 谐波平衡法是求解非线性动力学系统周期解的最常用方法, 但对非线性项进行高阶近似需要庞杂的公式推导, 限制了该方法的超高精度解算. 通过对频域非线性量的时域等价重构, 提出了重构谐波平衡法(RHB法), 解决了经典谐波平衡法超高阶次计算难题. 然而, 上述两种方法均要求动力学系统为多项式型非线性, 且无法直接用来求解非线性系统的拟周期解. 针对上述问题, 文章提出一种将RHB法和复杂非线性系统等价重铸法相结合的计算方法, 首先将一般非线性问题无损重铸为多项式型非线性系统, 然后用RHB法进行高精度求解; 针对拟周期响应求解问题, 提出基于“补频”思想的RHB方法, 通过基频的优化筛选, 实现拟周期响应的快速精准捕捉. 选取非线性单摆、相对论谐振子和非线性耦合非对称摆等典型系统进行仿真计算, 仿真结果表明, 所提出的RHB-重铸法在解非多项式型非线性系统的稳态响应时精度保持为${10^{ - 12}}$量级, 达计算机精度, 远超现有方法水平. 补频RHB法则实现了对拟周期问题的高效解算, 拓宽了方法对真实物理响应的求解范围.
-
关键词:
- 非多项式型非线性系统 /
- 重构谐波平衡法 /
- 微分方程重铸 /
- 非线性单摆 /
- 拟周期响应
Abstract: The harmonic balance method is the most commonly used method for solving periodic solutions of nonlinear dynamic systems, but the high-order approximation of nonlinear terms requires sophisticated formula derivations, which limits its ultra-high accuracy computation. The authors' team proposed the reconstruction harmonic balance (RHB) method through the equivalent reconstruction of the frequency domain nonlinear quantity in the time domain, which successfully conquered the problem of ultra-high-order calculation of the classical harmonic balance method. However, both two methods require the dynamical system to be polynomial nonlinear, and cannot be directly used to solve the quasi-periodic solution of the nonlinear system. In view of the above problems, this paper proposes a novel method that combines the RHB method and the recast technique for complex nonlinear systems. First, the general nonlinear problem is non-destructively recast into a polynomial nonlinear system, and then the RHB method is used for high-precision solutions. Aiming at computing the quasi-periodic response, which is hard to obtain relatively accurate results by considering only one base frequency, the RHB method based on the idea of "supplemental frequency" is derived. By optimizing and selecting base frequencies, the fast and accurate capture of quasi-periodic response is achieved in this paper. The typical systems such as the nonlinear pendulum, and the nonlinear coupling asymmetric pendulum are selected for simulation and algorithm verification. Simulation results show that the accuracy of the proposed RHB-recast method achieves an accuracy up to 10−12 when solving the non-polynomial type nonlinear systems, reaching the computer round-off-error accuracy, far exceeding the state-of-the-art methods. The supplemental frequency RHB method has been shown to produce efficient solutions for quasi-periodic problems, expanding the scope of the RHB-like methods for solving complex physical responses. -
引 言
工程中多数动力学问题, 其数学模型都是非线性的, 线性系统只是在一定假设及限制条件下对非线性系统的理想化近似[1]. 随着控制科学与航空宇航任务日益复杂, 强非线性的影响越发不容忽视. 谐波平衡法(harmonic balance method, HB)在求解非线性系统周期解中应用广泛, 但作为一种半解析半数值方法, 随着系统自由度、方法阶数的增加, 公式推导工作将变得困难[2]. 借助计算机数学软件可辅助推导工作, 但对计算效率及性能的提升作用仍然有限. 采用7阶HB法推导立方非线性项, Mathematica代码长达1321行. Hall等[3]提出了高维谐波平衡法(high-dimensional harmonic balance, HDHB), 通过“频域非线性量的时域近似计算”来简化公式推导, 但是, 由于引入近似关系导致非物理假解问题[4-5]. Dai等[6]发现了频域和时域非线性项之间的条件等价恒等式, 基于此, 首创了重构谐波平衡法(reconstruction harmonic balance, RHB)实现超高阶(N > 100)高精计算, 并给出时域配点计算的最优采样定理, 从理论上消除非物理解.
HB类方法(HB法及其改进方法)不仅用于计算简单非线性系统(杜芬、范德波方程等), 还发展到航空航天、深空探测等前沿领域. 哈尔滨工业大学陈毅等[7]使用谐波平衡−交变频域/时域法(HB-alternating frequency-time, HB-AFT)用于求解航空发动机中双转子−轴承−机匣系统动力学方程. 中山大学陈衍茂等[8]使用增量谐波平衡法对带机外挂载的二元机翼进行动力学特性研究. RHB法也被作为三体轨道的设计依据, 计算结果符合鹊桥中继卫星实际飞行数据[6].
但是该方法仍面临着两个问题: (1)谐波平衡法的本质依赖于非线性项的傅里叶级数展开, 因此, 受限于多项式型非线性系统求解, 对于非多项式型复杂非线性问题, 谐波平衡法难以适用; (2)已有谐波平衡类方法建立在单基频的假设上进行级数展开, 由于拟周期响应存在多基频的特征, 因此已有方法难以直接求解高精度拟周期解.
针对非多项式型非线性系统的求解问题, 目前有两类处理方法, 分别为直接法和间接法. 直接法包括HB-Taylor法[9]与HB-AFT法[10-11]. HB-Taylor法通过泰勒级数展开将非线性函数用有限阶近似多项式描述; HB-AFT法通过对非线性项的时域值采样以离散原问题. 由于直接法对原系统进行了近似, 导致求解精度低, 且计算性能分别受制于高阶级数描述和过采样等问题. Cochelin等[12]提出的谐波平衡−重铸法(HB-recast)是一种间接方法, 通过重铸(recast)技术, 成功将复杂非线性微分动力学系统无损变换为多项式型微分代数方程, 然后用HB法加以求解[13-14]. 但是, 受限于HB法的高阶计算(使用重铸法会增加系统的维度, 进一步增加了高阶计算的难度)与原重铸形式的二次型限制(原方法要求新系统中非线性至多为二次多项式), 至今难以对复杂非线性系统周期响应进行高效高精求解.
第2个难题是拟周期响应求解问题. 拟周期响应大量出现在非线性动力学系统中[15-17], 其频率响应由多个不可约基频及其线性组合描述[18]. 由于传统HB类方法基于单个基频进行近似解逼近, 不能简单通过基频及其整数倍频率分量对拟周期响应描述, Chua等[19]提出了结合广义傅里叶级数的改进多频HB法, 实现拟周期响应的准确求解. 然而, 使用多频HB类方法对强非线性项进行频域内高阶描述时, 计算效率严重受限于高维傅里叶分析(频域分量由多重积分[20]、求和[21]计算得到). Liu等[22]使用多频HDHB法求解受迫范德波振子的稳态响应, 避免非线性项谐波系数的直接表示以提高求解效率, 但是多频HDHB法的精度受到非物理解的破坏[23]. 总之, 当前关于HB类方法的研究已拓宽到对拟周期响应求解领域, 然而由于多基频计算中的高阶频域描述困难, 对此类复杂响应的高性能求解尚待解决.
针对上述问题, 本文提出了RHB法与重铸法、多频谐波平衡计算相结合的两种方法: (1) RHB-重铸法; (2)重构多谐波平衡法(reconstruction multiple harmonic balance, RMHB). 一方面RHB-重铸法通过将一般非线性等价转化为多项式型非线性系统, 再采用RHB法以确定时域最优采样点数, 实现对复杂非线性动力学系统的高阶高精求解, 仿真误差达计算机精度. 另一方面, RMHB法通过筛选和补充多个基频, 借鉴RHB法的时域等价重构思想, 完善了学术界在使用多频HB类方法求解拟周期响应时消除混淆假解的理论研究.
1. 谐波平衡法及重构谐波平衡法
1.1 谐波平衡法
对非线性动力学系统
$$ {\boldsymbol{\dot x}} = {\boldsymbol{f}}({\boldsymbol{x}},t) $$ (1) HB法用有限阶傅里叶级数来构造近似解及其导数
$$ {\boldsymbol{x}}\left( t \right) = {{\boldsymbol{x}}_0} + \sum\limits_{n = 1}^N {\left( {{{\boldsymbol{x}}_{2n - 1}}\cos n\omega t + {{\boldsymbol{x}}_{2n}}\sin n\omega t} \right)} $$ (2) $$ {\boldsymbol{\dot x}}\left( t \right) = \sum\limits_{n = 1}^N {\left( { - n\omega {{\boldsymbol{x}}_{2n - 1}}\sin n\omega t + n\omega {{\boldsymbol{x}}_{2n}}\cos n\omega t} \right)} $$ (3) 其中$N$是HB法的阶数. ${{\boldsymbol{x}}_0},{{\boldsymbol{x}}_1}, \cdots ,{{\boldsymbol{x}}_{2 n}}$为未知傅里叶系数, 也被称为频域变量, 将频域变量组成待求解向量${\boldsymbol{\hat x}} = {[{{\boldsymbol{x}}_0}{\text{ }}{{\boldsymbol{x}}_1}{\text{ }} \cdots {\text{ }}{{\boldsymbol{x}}_{2 n}}]^{\rm{T}}}$. 假设多项式函数$f(x) = {x^\phi },$ 忽略由非线性项而出现的高次谐波, 只需展开前$N$次的傅里叶分量
$$ f\left( {x,t} \right) = {x^\phi }\left( t \right) = {r_0} + \sum\limits_{n = 1}^N {\left( {{r_{2n - 1}}\cos n\omega t + {r_{2n}}\sin n\omega t} \right)} $$ (4) 其中各次分量${r_0},{r_1}, \cdots ,{r_{2 n}}$为
$$ \left. \begin{split} & {{r_0} = \frac{1}{{2\text{π} }}\int_0^{2\text{π} } {{{\left[ {{x_0} + \sum\limits_{n = 1}^N {\left( {{x_{2n - 1}}\cos n\theta + {x_{2n}}\sin n\theta } \right)} } \right]}^\phi }} {\rm{d}}\theta } \\ & {r_{2n - 1}} = \frac{1}{\text{π} }\int_0^{2\text{π} } \bigg[ {x_0} + \sum\limits_{n = 1}^N \left( {x_{2n - 1}}\cos n\theta +\right. \\ &\qquad \left.{x_{2n}}\sin n\theta \right) \bigg]^\phi \cos {n\theta } {\rm{d}}\theta \\ & {r_{2n}} = \frac{1}{\text{π} }\int_0^{2\text{π} } \bigg[ {x_0} + \sum\limits_{n = 1}^N \left( {x_{2n - 1}}\cos n\theta +\right.\\ &\qquad \left.{x_{2n}}\sin n\theta \right) \bigg]^\phi \sin {n\theta } {\rm{d}}\theta \end{split} \right\} $$ (5) 其中$n = 1,2, \cdots ,N$; $\theta = \omega t$. 将表示非线性函数的傅里叶分量记为${\boldsymbol{\hat f}}$, 是待求解向量${\boldsymbol{\hat x}}$的非线性函数.
将式(2) ~ 式(4)代入微分方程(1), 令常数项及前$n$次谐波$\cos n\omega t$和$\sin n\omega t\left( {n = 1,2, \cdots ,N} \right)$系数配平, 得到非线性代数方程组
$$ \omega {\boldsymbol{A\hat x}} = {\boldsymbol{\hat f}}({\boldsymbol{\hat x}}) $$ (6) 其中, 分块矩阵${\boldsymbol{A}}$为
$$ {\boldsymbol{A}} = {\text{diag}}\left( {\left[ {0,{{\boldsymbol{J}}_1},{{\boldsymbol{J}}_2}, \cdots ,{{\boldsymbol{J}}_N}} \right]} \right) \text{, } {{\boldsymbol{J}}_n} = n\left[ {\begin{array}{*{20}{c}} 0&1 \\ { - 1}&0 \end{array}} \right] $$ 由三角函数的求导关系得到.
非线性项的近似表示(4)中, 符号运算的复杂度会随算法阶次的提高呈指数增长[6]. 当HB法的阶数超过20, 即便有计算机数学软件的辅助, 非线性项的推导整理工作量仍难以接受.
1.2 高维谐波平衡法
为简化HB法的符号运算量, Hall等[3]使用时域值替代频域分量, 即将$N$阶HB法中的傅里叶系数作用转换矩阵${{\boldsymbol{E}}_{{\rm{HDHB}}}}$, 建立与一个周期$2 N + 1$个等距配点上时域量间的联系, 定义
$$ {\boldsymbol{\hat x}} = {{\boldsymbol{E}}_{{\text{HDHB}}}}{\boldsymbol{\tilde x}},\quad {\boldsymbol{\hat f}} \approx {{\boldsymbol{E}}_{{\text{HDHB}}}}{\boldsymbol{\tilde f}} $$ (7) 其中
$$ {\boldsymbol{\tilde x}} = {[{\boldsymbol{x}}\left( {{t_0}} \right){\text{ }}{\boldsymbol{x}}\left( {{t_1}} \right){\text{ }} \cdots {\text{ }}{\boldsymbol{x}}\left( {{t_{2N}}} \right)]^{\text{T}}} $$ $$ {\boldsymbol{\tilde f}} = {[{\boldsymbol{f}}\left( {{\boldsymbol{x}}\left( {{t_0}} \right)} \right){\text{ }}{\boldsymbol{f}}\left( {{\boldsymbol{x}}\left( {{t_1}} \right)} \right){\text{ }} \cdots {\text{ }}{\boldsymbol{f}}\left( {{\boldsymbol{x}}\left( {{t_{2N}}} \right)} \right)]^{\text{T}}} $$ $$ \begin{split} &{{\boldsymbol{E}}_{{\text{HDHB}}}} = \frac{1}{{2N + 1}}\\ &\quad\left[ {\begin{array}{*{20}{c}} {{1 / 2}}&{{1 / 2}}& \cdots &{{1 / 2}} \\ {\cos \omega {t_0}}&{\cos \omega {t_1}}& \cdots &{\cos \omega {t_{2N}}} \\ {\sin \omega {t_0}}&{\sin \omega {t_1}}& \cdots &{\sin \omega {t_{2N}}} \\ \vdots & \vdots &\vdots & \vdots\\ {\cos N\omega {t_0}} & {\cos N\omega {t_1}}&\cdots&{\cos N\omega {t_{2N}}}\\ {\sin N\omega {t_0}}& {\sin N\omega {t_1}}&\cdots&{\sin N\omega {t_{2N}}} \end{array}} \right] \end{split}$$ 将式(7)代入HB法代数方程组(6)得到HDHB法求解微分方程(1)对应的代数方程组
$$ \omega {\boldsymbol{A}}{{\boldsymbol{E}}_{{\text{HDHB}}}}{\boldsymbol{\tilde x}} = {{\boldsymbol{E}}_{{\text{HDHB}}}}{\boldsymbol{\tilde f}}({\boldsymbol{\tilde x}}) $$ (8) HDHB法显著提高计算效率, 并被认为是HB法的一种改进, 在计算流体力学领域得到应用[2,24]. 但该算法求解强非线性动力学问题时产生非物理解(也称假解), 严重影响求解精度. 如图1所示, HDHB法计算结果虽然是方程组(8)的数学解(使求解算法收敛), 但已偏离了真实的物理响应.
Liu等[4]指出, 假解现象产生于对非线性项近似的过程中, 存在高阶谐波对低次的混淆(污染). 以立方项非线性${x^3}$为例, 二阶HB法中${\boldsymbol{\hat f}}({\boldsymbol{\hat x}})$表达式
$$ {{\boldsymbol{\hat f}}_{{\text{HB}}}} = {\left[ {\begin{array}{*{20}{c}} {\begin{array}{*{20}{c}} {{r_0}}&{{r_1}} \end{array}}&{{r_2}}&{{r_3}}&{{r_4}} \end{array}} \right]^{\text{T}}} $$ HDHB法计算中, 非线性项的频域分量(9)中不仅包含了原HB法中的所有项, 还包括附加杂项, 这些杂项的混入导致求解中出现非物理解
$$ {{\boldsymbol{\hat f}}_{{\rm{HDHB}}}} = \left[ {\begin{array}{*{20}{c}} {{r_0} + \dfrac{3}{4}\left( {x_3^2 - x_4^2} \right){x_1} - \dfrac{3}{2}{x_2}{x_3}{x_4}} \\ {{r_1} + \dfrac{3}{4}\left( {x_1^2 - x_2^2} \right){x_3} + \dfrac{3}{2}\left( {x_3^2 - x_4^2} \right){x_0} + \dfrac{1}{4}x_3^3 - \dfrac{3}{4}{x_3}x_4^2 - \dfrac{3}{2}{x_1}{x_2}{x_4}} \\ {{r_2} - \dfrac{3}{4}\left( {x_1^2 - x_2^2} \right){x_4} - \dfrac{1}{4}x_4^3 + \dfrac{3}{4}x_3^2{x_4} - \dfrac{3}{2}{x_1}{x_2}{x_3} - 3{x_0}{x_3}{x_4}} \\ {{r_3} + \dfrac{1}{4}x_1^3 + \dfrac{3}{4}\left( {x_3^2 - x_4^2} \right){x_1} - \dfrac{3}{4}{x_1}x_2^2 + \dfrac{3}{2}{x_2}{x_3}{x_4} + 3{x_0}{x_1}{x_3} - 3{x_0}{x_2}{x_4}} \\ {{r_4} + \dfrac{1}{4}x_2^3 + \dfrac{3}{4}\left( {x_3^2 - x_4^2} \right){x_2} - \dfrac{3}{4}x_1^2{x_2} - \dfrac{3}{2}{x_1}{x_3}{x_4} - 3{x_0}{x_2}{x_3} - 3{x_0}{x_1}{x_4}} \end{array}} \right] $$ (9) 1.3 重构谐波平衡法
Dai等[25]证明, HDHB法与时域配点法等价, 时域配点法通过使用时域配点上函数值对微分方程进行离散, 从而联立代数方程组
$$ \omega {\boldsymbol{EA\hat x}} = {\boldsymbol{\tilde f}}({\boldsymbol{\tilde x}}) $$ (10) 其中配点矩阵
$$ \begin{split} &{\boldsymbol{E}} = \\ &\left[ {\begin{array}{*{20}{c}} 1&{\cos \omega {t_0}}&{\sin \omega {t_0}}& \cdots &{\cos N\omega {t_0}}&{\sin N\omega {t_0}} \\ 1&{\cos \omega {t_1}}&{\sin \omega {t_1}}& \cdots &{\cos N\omega {t_1}}&{\sin N\omega {t_1}} \\ \vdots & \vdots & \vdots &{\vdots}& \vdots & \vdots \\ 1&{\cos \omega {t_M}}&{\sin \omega {t_M}}& \cdots &{\cos N\omega {t_M}}&{\sin N\omega {t_M}} \end{array}} \right] \end{split}$$ 混淆是造成计算出现非物理解的原因. 至于高次谐波在时域配点计算中影响低次谐波的机理, 则遵循“混淆规则”[25].
混淆规则: 假设$\alpha \in \left[ {0,\text{π} } \right]$被以$h$为间隔均分为离散时间点, 那么配点法中可区分的最大谐波次数$n \in \left[ { - L,L} \right]$, 其中$L = {\text{π} \mathord{\left/ {\vphantom {\text{π} h}} \right. } h}$, $L$为“极限波次”. 高次谐波$n\left( {\left| n \right| > L} \right)$被误认为是相应低次谐波${n_a}$
$$ {n_a} = n - 2mL $$ 其中${n_a} \in \left[ { - L,L} \right]$, $m$是整数.
混淆规则指出, 时域配点法可区分的最大谐波次数由配点数(离散时间点的间隔)决定. 增加时域配点法中的配点数量, 可区分的谐波次数越高, 此时方程(10)个数多于待求解变量数, 配点矩阵${\boldsymbol{E}}$为列满秩矩阵, 存在Moore-Penrose广义逆矩阵$ {{\boldsymbol{E}}^*} $
$$ {{\boldsymbol{E}}^*} = \frac{2}{M}{\left[ {\begin{array}{*{20}{c}} {{1 \mathord{\left/ {\vphantom {1 2}} \right. } 2}}&{{1 \mathord{\left/ {\vphantom {1 2}} \right. } 2}}& \cdots &{{1 \mathord{\left/ {\vphantom {1 2}} \right. } 2}} \\ {\cos \omega {t_1}}&{\cos \omega {t_2}}& \cdots &{\cos \omega {t_M}} \\ {\sin \omega {t_1}}&{\sin \omega {t_2}}& \cdots &{\sin \omega {t_M}} \\ \vdots & \vdots &{\vdots}& \vdots \\ {\sin N\omega {t_1}}&{\sin N\omega {t_2}}& \cdots &{\sin N\omega {t_M}} \end{array}} \right]_{\left( {2N + 1} \right) \times M}} $$ 使得
$$ {{\boldsymbol{E}}^*}{\boldsymbol{E}} = {{\boldsymbol{I}}_{2N + 1}} $$ 其中$ {{\boldsymbol{I}}_{2 N + 1}} $为单位矩阵. 在配点法方程同时作用矩阵$ {{\boldsymbol{E}}^*} $, 对方程组降维, 得到RHB法代数方程组
$$ \omega {\boldsymbol{A\hat x}} = {{\boldsymbol{E}}^*}{\boldsymbol{\tilde f}}({\boldsymbol{\tilde x}}) $$ (11) RHB法在保证算法效率的同时, 消除非物理解, 从而实现对原HB法的最佳重构. 定理1围绕如何选择合适的配点数, 给出消除混淆确定条件[6].
定理1(条件等价性): 如果配点数$M$, 方法阶数$N$, 多项式非线性的幂次$\phi $满足
$$ M > \left( {\phi + 1} \right)N $$ (12) 则RHB法与HB法等价.
为说明RHB法消除混淆的效果. 分别使用HDHB法和RHB法计算杜芬方程, 在分岔处(频率$\omega = 2$)进行蒙特卡洛法模拟, 选取104组随机初值(各频域分量在区间[−2, 2]中选取), 得到如表1所示的统计结果. 统计结果表明, HDHB法计算产生58组解, 其中55组都是非物理解; 而RHB法只计算得到3个具有物理含义的解[6].
表 1 RHB法与HDHB法在分岔处的解分布Table 1. Statistical distribution of solution by the RHB and the HDHB method at bifurcation pointMethod Upper
branch/%Lower
branch/%Unstable
branch/%Non-physical/% RHB 57.13 19.71 23.16 0 HDHB 29.14 14.21 15.96 48.69 此外, HB-AFT法与HDHB法计算流程一致, 但HB-AFT法的原理是, 选用配点数$2\phi N + 1$来消除混淆[11], 但过采样会占用计算机资源, 在实际计算时会导致更大的CPU和RAM计算负担[2]. RHB法基于时域配点法的统一框架, 根据配点数的差异将所有HB类方法(HDHB法、HB-AFT法等)建立起联系.
2. 重构谐波平衡法改进策略
2.1 微分方程重铸技术
针对非多项式型非线性系统的HB法求解, Cochelin等[12]提出将原系统改写为二次型系统
$$ {\boldsymbol{\dot z}} = {\boldsymbol{c}} + {\boldsymbol{l}}\left( {\boldsymbol{z}} \right) + {\boldsymbol{q}}\left( {{\boldsymbol{z}},{\boldsymbol{z}}} \right) $$ (13) 其中未知向量${\boldsymbol{z}}$包含微分方程的原始变量${\boldsymbol{x}}$及一些新的变量(引入的新变量用以改写方程). ${\boldsymbol{c}}$是关于未知量${\boldsymbol{z}}$的常数向量; ${\boldsymbol{l}}$是关于${\boldsymbol{z}}$的线性向量值运算符; ${\boldsymbol{q}}$则是二次向量值运算符.
因为对任意${x^\phi }$次非线性, RHB法都能实现高效计算, 改写后方程可以是更高次多项式, 将式(13)进一步写成适配于RHB法的微分方程重铸形式
$$ {\boldsymbol{\dot z}} = {\boldsymbol{c}} + {\boldsymbol{l}}\left( {\boldsymbol{z}} \right) + {\boldsymbol{n}}\left( {\boldsymbol{z}} \right) $$ (14) 其中${\boldsymbol{n}}$可以是关于变量${\boldsymbol{z}}$的任意$\phi $次多项式函数. 重铸格式(14)涵盖多类型非多项式型非线性问题, 下面将分类加以介绍.
(1)微分项耦合型
对范德波方程分析, 一阶导数$\dot x$与平方项${x^2}$耦合
$$ \ddot x - \varepsilon \left( {1 - {x^2}} \right)\dot x + x = F\cos \omega t $$ (15) 通过将方程重铸为标准形式(14), 得到
$$ \begin{split} & {\dot x = u} \\ & {\dot u = \left( {\varepsilon u - x} \right) - \varepsilon {x^2}u + F\cos \omega t} \end{split} $$ 原方程转化为典型的非线性度$\phi = 3$的动力学系统. 得出结论: 导数项的耦合不影响非线性度计入, 在实际HB计算中, 任意阶导数只计入一个非线性度.
(2)有理分式型
以Rayleigh-Plesset方程为例
$$ R\ddot R = - \frac{3}{2}{\dot R^2} - A\frac{{\dot R}}{R} - B\frac{1}{R} + C\frac{1}{{{R^3}}} + D - E\cos {\omega t} $$ (16) 式中, $A$, $B$, $C$, $D$均为常系数. 将方程两端同除以$R$, 并运用倒数关系, 令$v = \dot R$, $x = {1 \mathord{\left/ {\vphantom {1 R}} \right. } R}$得到[6,12]
$$\begin{split} & {\dot R = v} \\ & {\dot v = - \frac{3}{2}{v^2}x - Av{x^2} - B{x^2} + C{x^4} + Dx - Ex\cos {\omega t} } \\ & {0 = Rx - 1} \end{split} $$ 处理有理分式型非线性项, 通过额外增加方程$Rx = 1$, 将方程改写为非线性度$\phi = 4$的多项式型非线性系统.
(3)根式型
相对论谐振子方程(17)中非线性项为根式
$$ \ddot x + {\left( {1 - {{\dot x}^2}} \right)^{{3 \mathord{\left/ {\vphantom {3 2}} \right. } 2}}}x = 0 $$ (17) 引入新变量${u^2} = 1 - {\dot x^2}$改写根号内关系式, 利用根式的特性${\left( {1 - {{\dot x}^2}} \right)^{{3 \mathord{\left/ {\vphantom {3 2}} \right. } 2}}} = {\left( {{u^2}} \right)^{{3 \mathord{\left/ {\vphantom {3 2}} \right. } 2}}} = {u^3}$, 将方程改写为非线性度$\phi = 4$的多项式型非线性系统
$$ \left. \begin{split} & {\dot x - v = 0} \\ & {\dot v + {u^3}x = 0} \\ & {{u^2} + {v^2} - 1 = 0} \end{split}\right\} $$ (18) 对任意形如${\alpha ^{{q \mathord{\left/ {\vphantom {q p}} \right. } p}}}$的根式, 令${\beta ^p} = \alpha $, 使原根式型非线性转化得到多项式型: ${\alpha ^{{q \mathord{\left/ {\vphantom {q p}} \right. } p}}} = {\left( {{\beta ^p}} \right)^{{q \mathord{\left/ {\vphantom {q p}} \right. } p}}} = {\beta ^q}$.
(4)初等超越函数
对初等超越函数的处理, 需要导数信息来实现对原微分方程的重铸. 以非线性单摆方程为例
$$ \ddot \theta \left( t \right) + \sin \theta \left( t \right) = 0 $$ (19) 由于三角函数是超越函数, 不能通过简单的代数关系式改写转化为多项式型. 首先额外引入两个变量$s$, $c$来分别表示$ s\left( t \right) = \sin \theta \left( t \right) $, $c\left( t \right) = \cos \theta \left( t \right)$. 利用三角函数的导数性质
$$ \begin{gathered} \frac{{\text{d}}}{{{\text{d}}t}}\sin \theta \left( t \right) = \cos \theta \left( t \right) \cdot \dot \theta \left( t \right) \\ \frac{{\text{d}}}{{{\text{d}}t}}\cos \theta \left( t \right) = - \sin \theta \left( t \right) \cdot \dot \theta \left( t \right) \\ \end{gathered} $$ 建立补充方程, 从而实现对微分方程组的重铸[13]
$$ \left. \begin{split} & {\dot \theta \left( t \right) = v\left( t \right)} \\ & {\dot v\left( t \right) = - s\left( t \right)} \\ & {\dot s\left( t \right) = c\left( t \right)v\left( t \right)} \\ & {\dot c\left( t \right) = - s\left( t \right)v\left( t \right)} \end{split} \right\} $$ (20) 初等超越函数非线性的重铸, 以导数关系作为方程组改写的依据(部分需要用到有理分式、根式型非线性重铸技巧). 附录表格中罗列了几类常见初等超越函数改写思路与重铸形式, 可供参考.
总之, 本文主要针对非线性项为连续函数情形(包括微分项耦合、有理分式型、根式型以及初等超越函数4类), 重铸技术通过变量替换、利用特殊函数的导数关系(见附表1)等关键步骤将一般连续函数非线性项改写为更利于重构谐波平衡法求解的多项式形式(14).
2.2 多谐波平衡计算
“补频”(supplemental frequency)思想[17]是在原来单频假设的基础上, 补充多个与响应相关的频率. 假设考虑两个基频, 将假设函数改写为
$$ \begin{split} & {x_i}\left( t \right) = \sum\nolimits_m^{} {\sum\nolimits_n^{} {{x_{ic}}\left( {m,n} \right)\cos \left[ {\left( {m{\omega _1} + n{\omega _2}} \right)t} \right]} }+ \\ &\qquad {x_{is}}\left( {m,n} \right)\sin \left[ {\left( {m{\omega _1} + n{\omega _2}} \right)t} \right] \end{split} $$ (21) 参数$m$和$n$满足不等式
$$ \left|m\right| + \left|n\right|\leqslant p $$ 其中$p$代表二维傅里叶级数的截断[19-20], 类似于RHB法中阶数$N$.
以RHB法理论为基础, 提出的RMHB法利用时域信息, 将动力学问题(1)转化为非线性代数方程(11)求解. 多基频计算中, 配点矩阵${\boldsymbol{E}}$以及转换矩阵${{\boldsymbol{E}}^*}$形式, 记${c^{m,n}} = \cos (m{\omega _1} + n{\omega _2})t$, ${s^{m,n}} = \sin (m{\omega _1} + n{\omega _2})t$. 在RMHB法计算中, 存在选取合适的采样周期和时域配点数用以消除混淆的结论[23].
定理2(多频计算中条件等价性): 假设多项式非线性的幂次$\phi $的非线性系统响应中包含两个基频, 基频之比${{{\omega _1}} \mathord{\left/ {\vphantom {{{\omega _1}} {{\omega _2}}}} \right. } {{\omega _2}}}$为有理数. 则RMHB法消除混淆需满足采样时间$ T = {{2\text{π} } \mathord{\left/ {\vphantom {{2\text{π} } {{\text{GCD}}\left( {{\omega _1},{\omega _2}} \right)}}} \right. } {{\text{GCD}}\left( {{\omega _1},{\omega _2}} \right)}} $, 配点数
$$ M > \left( {\phi + 1} \right)\frac{{p \cdot \max \left( {{\omega _1},{\omega _2}} \right)}}{{{\text{GCD}}\left( {{\omega _1},{\omega _2}} \right)}} $$ (22) 其中${\rm{GCD}}$为两数最大公因数.
$$ \left.\begin{split} & {{\boldsymbol{E}} = {{\left[ {\begin{array}{*{20}{c}} 1&1& \cdots &1 \\ {{c^{1,0}}\left( {{t_1}} \right)}&{{c^{1,0}}\left( {{t_2}} \right)}& \cdots &{{c^{1,0}}\left( {{t_M}} \right)} \\ {{s^{1,0}}\left( {{t_1}} \right)}&{{s^{1,0}}\left( {{t_2}} \right)}& \cdots &{{s^{1,0}}\left( {{t_M}} \right)} \\ \vdots & \vdots & \vdots & \vdots \\ {{c^{m,n}}\left( {{t_1}} \right)}&{{c^{m,n}}\left( {{t_2}} \right)}& \cdots &{{c^{m,n}}\left( {{t_M}} \right)} \\ {{s^{m,n}}\left( {{t_1}} \right)}&{{s^{m,n}}\left( {{t_2}} \right)}& \cdots &{{s^{m,n}}\left( {{t_M}} \right)} \\ \vdots & \vdots & \vdots & \vdots \\ {{c^{0,p}}\left( {{t_1}} \right)}&{{c^{0,p}}\left( {{t_2}} \right)}& \cdots &{{c^{0,p}}\left( {{t_M}} \right)} \\ {{s^{0,p}}\left( {{t_1}} \right)}&{{s^{0,p}}\left( {{t_2}} \right)}& \cdots &{{s^{0,p}}\left( {{t_M}} \right)} \end{array}} \right]}^{\text{T}}}} \\ & {{{\boldsymbol{E}}^*} = \frac{2}{M}\left[ {\begin{array}{*{20}{c}} {\dfrac{1}{2}}&{\dfrac{1}{2}}& \cdots &{\dfrac{1}{2}} \\ {{c^{1,0}}\left( {{t_1}} \right)}&{{c^{1,0}}\left( {{t_2}} \right)}& \cdots &{{c^{1,0}}\left( {{t_M}} \right)} \\ {{s^{1,0}}\left( {{t_1}} \right)}&{{s^{1,0}}\left( {{t_2}} \right)}& \cdots &{{s^{1,0}}\left( {{t_M}} \right)} \\ \vdots & \vdots & \vdots & \vdots \\ {{c^{m,n}}\left( {{t_1}} \right)}&{{c^{m,n}}\left( {{t_2}} \right)}& \cdots &{{c^{m,n}}\left( {{t_M}} \right)} \\ {{s^{m,n}}\left( {{t_1}} \right)}&{{s^{m,n}}\left( {{t_2}} \right)}& \cdots &{{s^{m,n}}\left( {{t_M}} \right)} \\ \vdots & \vdots & \vdots & \vdots \\ {{c^{0,p}}\left( {{t_1}} \right)}&{{c^{0,p}}\left( {{t_2}} \right)}& \cdots &{{c^{0,p}}\left( {{t_M}} \right)} \\ {{s^{0,p}}\left( {{t_1}} \right)}&{{s^{0,p}}\left( {{t_2}} \right)}& \cdots &{{s^{0,p}}\left( {{t_M}} \right)} \end{array}} \right]} \end{split} \right\} $$ (23) 3. 数值仿真结果
3.1 相对论谐振子
作为物理学中经典问题, 有必要对谐振子模型进行完整而严格的相对论处理[26]. 考虑一个静质量为$m$的质点在一维谐振力$F = - k\bar x$的作用下做相对论运动. 其中$k$为弹性常数, $\bar x$为位移量. 结合牛顿运动运动学方程以及动量定理, 可以推导得到相对论振子方程[27-28]
$$ \frac{{{{\text{d}}^2}\bar x}}{{{\text{d}}{{\bar t}^2}}} + \frac{k}{m}{\left[ {1 - \frac{1}{{{c^2}}}{{\left( {\frac{{{\text{d}}\bar x}}{{{\text{d}}\bar t}}} \right)}^2}} \right]^{{3 \mathord{\left/ {\vphantom {3 2}} \right. } 2}}}x = 0 $$ 其中$\bar t$是时间坐标(维度变量), 对方程进行无量纲化得到方程(17). 初值条件$x\left( 0 \right) = 0$, $\dot {x}\left( 0 \right) = \beta $, 其中$ - 1 < \beta < 1$, 且非线性振子的振动周期与对应周期解依赖于初始速度$\beta $. Mickens为了使相轨迹充满整个相平面引入非线性变换[29]
$$ y = \frac{w}{{\sqrt {1 + {w^2}} }} $$ 再使用HB法计算. 该处理方式在谐波振子低速运动时提供了高精度解[29-30]. 当振子高速运动接近光速时($\beta = 0.85$), 如图2所示, 按照Mickens变换法求解将产生较大计算偏差.
RHB-重铸法按照根式型非线性重铸原则, 将相对论谐振子改写为方程(18), 再使用高阶RHB法求解. 结合表2和图3, RHB-重铸法可对高速运动的谐振子($\beta = 0.85$)直接进行高阶计算. 55阶RHB-重铸法能将计算误差控制在${10^{ - 12}}$量级(以数值积分为参照, 全文采用MATLAB内置函数ode15 s, 相对精度容差${10^{ - 13}}$, 绝对误差容限${10^{ - 20}})$, 总计算耗时在1 s内.
表 2 各阶RHB-重铸法求解相对论谐振子Table 2. Solving relativistic harmonic oscillator by the RHB-recast method with different ordersOrder of method Amplitude error Computing time/s 25 $4.36 \times {10^{ - 7}}$ 0.73 35 $1.95 \times {10^{ - 9}}$ 0.77 55 $3.87 \times {10^{ - 12}}$ 1.07 除通过方法阶数对求解精度的提高, 还能通过两种方式实现对计算性能的改善.
(1) 非线性代数方程组求解算法
非线性代数方程(11)求解算法能提高HB类方法计算性能, Zheng等[31]结合Tikhonov正则化求解, 黄建亮等[32]引入狗腿法思想结合回溯线搜索算法求解, Thomas等[33]使用Broyden法以提高HB法计算性能. 本文根据算例来说明L-M (Levenberg-Marquardt)法较之传统迭代方法在求解RHB方程组时体现出的优势. 图4展示当$\beta = 0.85$, 频域量初值估计${x_2} = 1$, ${u_0} = 0.7$, 使用两种不同的方程组求解算法: 牛顿迭代 (Newton-Raphson)法与L-M法得到的一个周期内误差曲线. 不同于牛顿迭代法, L-M法迭代格式为[34]
$$ {{\boldsymbol{x}}_{k + 1}} = {{\boldsymbol{x}}_k} - {\left[ {{\boldsymbol{J}}_k^{\text{T}}{{\boldsymbol{J}}_k} + {\lambda _k}{\text{diag}}\left( {{\boldsymbol{J}}_k^{\text{T}}{{\boldsymbol{J}}_k}} \right)} \right]^{ - 1}}{\boldsymbol{J}}_k^{\text{T}}{{\boldsymbol{F}}_k} $$ 其中${{\boldsymbol{x}}_k}$为当前迭代未知变量值, ${{\boldsymbol{F}}_k}$为函数值, ${{\boldsymbol{J}}_k}$为雅可比矩阵. L-M法通过衡量每步迭代的误差是否发散, 来决定撤回迭代并将参数$ \lambda_{k} $按十倍放缩. L-M法在求解非线性方程组时显示出更高的计算精度, 同比牛顿迭代法, 精度提高了${10^4}$以上.
(2) 合理选择代数方程
相对论谐振子(17)为保守系统, 振动频率$\omega $是状态变量[35]. 使用RHB-重构法对$n$自由度保守系统求解时, 共需考虑$ n\left( {2 N + 1} \right) + 1 $个变量, 需额外增加一个代数方程使RHB方程组适定. 本文以初值约束来探讨方程对计算性能的影响. 图5为采用55阶RHB-重构法与L-M法求解器, 使用不同代数方程(约束条件)得到的误差曲线. 分别对应: 初始位移约束$x\left( 0 \right) = 0$; 初始速度约束$\dot x\left( 0 \right) = \beta $; 同时考虑初始位移与速度约束. 同时考虑位移与速度约束时, RHB-重构法计算精度更高.
相对论谐振子在高速运动时显示出复杂的动力学特性, 需要高阶($ > 20$阶)方法来进行分析. 使用重铸技术, 将根式非线性转化为多项式型, 再使用RHB法实现快速解算, 克服高阶估计的限制. 如图6所示, RHB-重铸法适用于不同初速度条件下相对论振子的计算.
3.2 非线性单摆
非线性单摆(19)是物理学中经典问题, 且时域响应具有典型的周期性. 但现有HB类方法不能做到高性能求解: HB-Taylor法需要采用高的截断阶来保证计算精度; HB-AFT法则需进行过采样以抑制混淆误差(如图7所示).
对非线性单摆重铸, 得到方程(20). 使用25阶RHB-重铸法求解(大角度摆动, $\theta \left( 0 \right) = 1.5$), 初值估计${x_1} = 1.423$, ${s_1} = 1.065$, ${c_0} = 1.028$, 辅助L-M法求解器, 计算结果如图8所示, 与数值积分参考解比较, 误差控制在${10^{ - 12}}$数量级以下.
分别采用牛顿迭代法、Tikhonov正则法与L-M法得到图9计算结果. 牛顿迭代法中雅可比阵奇异, 求解失效; Tikhonov法与L-M法计算所得的误差数量级相近, 但是由于Tikhonov法对正则参数的优化使求解更耗时[31]. 本例中, 使用L-M法计算RHB法方程组仅需0.72 s, 而Tikhonov法耗时达1.17 s, L-M法较之Tikhonov法效率提高了62.5%, 达到了相近的求解精度. 对比传统方法(牛顿迭代法)与优化方法(Tikhonov法), L-M法都更适合于RHB法方程组的求解.
非线性单摆问题是保守系统, 其周期解$\theta \left( t \right)$与响应频率都依赖于初始振幅$\theta \left( 0 \right)$. 使用RHB法求解重铸方程组(20)时, 本文给出3种代数方程组合方案, 以助于对比分析.
方案1: 对未知量$\theta $, $v$, $s$和$c$全谐波平衡(从常数项到N次谐波项), 计$ 4\left( {2 N + 1} \right) $多个方程. 增加初始振幅约束$\theta \left( 0 \right) = \alpha $.
方案2: 对未知量$\theta $和$v$全谐波平衡, 计$ 2\left( {2 N + 1} \right) $多个方程. $s$和$c$从1次谐波开始, 平衡系数到$N$次, 计$2 N$多个方程. 增加初始振幅约束$\theta \left( 0 \right) = \alpha $与两个对附加变量$s$, $c$的初值约束
$$ \left. \begin{split} & {s\left( 0 \right) = \sin \theta \left( 0 \right) = \sin \alpha } \\ & {c\left( 0 \right) = \cos \theta \left( 0 \right) = \cos \alpha } \end{split} \right\} $$ (24) 方案3: 对未知量$\theta $和$v$全谐波展开, 计$ 2\left( {2 N + 1} \right) $多个方程. $s$和$c$从1次谐波展开到$N$次, 计$2 N$多个方程. 增加初始速度约束$ v(0) = 0 $与两个对附加变量$s$, $c$的初值约束式(24).
分别采用3种方案得到的计算结果如图10所示. 对比方案2和方案3, 方案1的计算误差大, 主要原因是没有对附加变量$s$和$c$的初值进行约束. 方案3同时利用了初始位移与速度信息, 比方案2计算精度更高.
对比求解非线性单摆的3种方法(HB-AFT法、RHB-Taylor法和RHB-重铸法)计算结果. 考虑同阶截断$N = 25$, 综合图11误差曲线和表3的计算结果, RHB-重铸法确定最优时域配点数$M = 76$, 降低采样成本. RHB-Taylor法虽可采用高阶截断(15次多项式)保证计算精度, 但超越函数的级数表示降低计算效率, 计算用时达1.80 s. RHB-重铸法与HB-AFT法计算效率相当, 但将计算精度提高了两个数量级以上.
表 3 3种方法计算非线性单摆结果对比Table 3. Comparison of corresponding results of three methods for solving nonlinear pendulumMethods $M$ Average error Computing time/s RHB-recast 76 1.9 × 10−12 0.70 RHB-Taylor 401 8.9 × 10−10 1.82 HB-AFT 210 1.5 × 10−9 0.69 RHB-重铸法适用于求解初等超越函数非线性问题. 尤其是当单摆做大幅度振动时, 传统的线性化手段无法很好地捕捉动力学响应, 而RHB-重铸法则可以提供高精度解析解(图12).
3.3 非线性耦合的非对称摆
用绳将质点小球悬挂于一根固定在铁架上的弹性杆末端, 而弹性杆能在水平面与质点小球同步振动. 此物理模型在许多其他二维摆系统中普遍存在, 如球面摆、傅科摆等. 特别地, 当傅科摆由于不完美的悬挂或由于侧向运动引起非线性耦合而导致不对称, 可能会导致额外的旋转, 从而掩盖地面效应. 非线性耦合的非对称单摆微分形式可以写作[36]
$$ \left. \begin{split} & {\ddot x + \left( {1 - \kappa } \right)x = - \left( {1 - \kappa } \right)\left( {x{{\dot x}^2} + x{{\dot y}^2} + xy\ddot y + {x^2}\ddot x} \right)} \\ & {\ddot y + y = - \left( {y{{\dot y}^2} + y{{\dot x}^2} + xy\ddot x + {y^2}\ddot y} \right)} \end{split} \right\} $$ (25) 此方程中2阶导数耦合入非线性项, 导致使用传统方法进行直接求解变得棘手. 以传统的有限差分类方法(欧拉法、龙格库塔法、MATLAB内置ode算法等)而言, 需额外计算代数方程组来辅助求解. 而部分解析求解法的计算效果同样不佳, Jia等[36]曾采用多时间尺度法推导方程的解为
$$\left. \begin{split} & {x = a\cos {\sigma _x}t + \frac{{\left( {1 - \kappa } \right)a{b^2}}}{{2{b^2} - {a^2} - 4\kappa }}\cos \left( {2{\sigma _y} - {\sigma _x}} \right)t} \\ & {y = b\cos {\sigma _y}t + \frac{{\left( {1 - \kappa } \right){a^2}b}}{{2{a^2} - {b^2} + 4\kappa }}\cos \left( {2{\sigma _x} - {\sigma _y}} \right)t} \end{split} \right\}$$ (26) 令参数$\kappa = 0.01$, 初始条件$ x\left( 0 \right) = 0.1 $, $ y\left( 0 \right) = 0.2 $, $\dot x\left( 0 \right) = \dot y\left( 0 \right) = 0$.
以修正Chebyshev-Picard迭代 (modified Chebyshev-Picard iteration, MCPI)法计算结果为标准解[37-38]. 仿真得到如图13所示, 解(26)不仅推导复杂, 且与真实的物理解间有较大误差.
对多时间尺度解(26)分析得知, 动力学响应包含两个基频(所有频率成分都是基频的线性组合), 通过快速傅里叶变换(FFT)得到: ${\omega _1} = 0.985\;7$, ${\omega _2} = 0.993\;5$. 又从微分方程重铸的角度分析非对称摆微分方程(25), 因为高阶导数只算作一个非线性度, 该方程是具有3次多项式非线性的动力学系统. 使用5阶RMHB法进行计算得到图14所示仿真结果, 计算用时8.24 s, $x$方向振动幅值误差$9.21 \times {10^{ - 3}}$, $y$方向振动幅值误差$3.90 \times {10^{ - 7}}$.
对2阶导耦合型非线性系统, 有限差分方法不能直接求解, 一些解析求解方法的计算精度低. 考虑两个基频的RMHB法不仅保证计算效率, 还能对拟周期运动进行准确的预测.
4. 结 论
围绕复杂非线性动力学系统高效高精度周期解算需求, 本文基于微分方程重铸、“补频”等思想, 分别提出RHB-重铸法与RMHB法, 主要的工作与结论总结如下.
(1)提出RHB-重铸法, 将一般非线性问题改写为多项式型非线性方程, 配合RHB法确定消除非物理解所需的最优配点阈值, 实现对非多项式型非线性动力学系统周期响应进行高阶预测. 数值实验说明, 对比HB-AFT法的时域过采样和HB-Taylor法对非线性函数进行泰勒级数截断两种处理方式, RHB-重铸法的计算精度高达${10^{ - 12}}$, 且计算时间不超过1 s, 同时保证了求解的高精度与高效率.
(2)提出结合“补频”思想的RMHB法, 在原来单频假设的基础上, 补充多个与物理响应相关的频率. 借鉴RHB法中时域重构等价性, 给出多频谐波平衡计算的最优配点数结论, 并充分利用非线性量的时域采样代替复杂的高维傅里叶分析, 最终实现对拟周期响应的准确捕捉.
(3)通过解算相对论谐振子、非线性单摆以及非线性耦合的非对称摆问题验证了本文算法的有效性. 针对传统有限差分类方法求解困难的多自由度耦合系统, RHB及本文提出的两种方法都不受方程形式的限制. 此外, 在相同的系统参数设置下, 合理地选择代数方程求解器、代数约束方案, 可有效提高非线性问题的求解精度.
本文提出的RHB-重铸法与RMHB法对复杂非线性动力学系统周期性及拟周期响应的计算效率与精度相较于现行的方法与处理方式都具有优势. 在后续的研究中, 将探究RHB法在高精周期响应求解方面的工程应用.
数据可用性声明
支撑本研究的科学数据已在中国科学院科学数据银行 (science data bank) ScienceDB平台公开发布, 访问地址为: https://cstr.cn/31253.11.sciencedb.j00140.00022或https://doi.org/10.57760/sciencedb.j00140.00022.
附录A
A1 常见初等超越函数的重铸A1. Recast form of the most common elementary transcendental functionsOriginal function Differential relationship Companion variable Recast equation $u\left( t \right) = \exp {x\left( t \right)}$ $\dot u\left( t \right) = \dot x\left( t \right) \cdot u\left( t \right)$ $\dot u\left( t \right) = \dot x\left( t \right) \cdot u\left( t \right)$ $u\left( t \right) = {\log _a} {x\left( t \right)} $ $\dot u\left( t \right) = {{\dot x\left( t \right)} \mathord{\left/ {\vphantom {{\dot x\left( t \right)} {\left( {x\left( t \right) \cdot \ln a} \right)}}} \right. } {\left( {x\left( t \right) \cdot \ln a} \right)}}$ $ v\left( t \right) = {1 \mathord{\left/ {\vphantom {1 {x\left( t \right)}}} \right. } {x\left( t \right)}} $ $\left. \begin{aligned}& {\dot u\left( t \right) = {{\left( {v\left( t \right) \cdot \dot x\left( t \right)} \right)} / {\ln a}}} \\ & {x\left( t \right)v\left( t \right) - 1 = 0} \end{aligned}\right\}$ $\left. \begin{aligned}& {u\left( t \right) = \sin {x\left( t \right)} } \\ & {v\left( t \right) = \cos {x\left( t \right)} } \end{aligned} \right\}$ $\left.\begin{aligned}& {\dot u\left( t \right) = \dot x\left( t \right) \cdot \cos {x\left( t \right)}} \\ & {\dot v\left( t \right) = - \dot x\left( t \right) \cdot \sin {x\left( t \right)} } \end{aligned} \right\}$ $\left. \begin{aligned}& {\dot u\left( t \right) = \dot x\left( t \right) \cdot v\left( t \right)} \\ & {\dot v\left( t \right) = - \dot x\left( t \right) \cdot u\left( t \right)} \end{aligned}\right\}$ $u\left( t \right) = \tan {x\left( t \right)} $ $\dot u\left( t \right) = \dot x\left( t \right) \cdot \left( {1 + {u^2}\left( t \right)} \right)$ $\dot u\left( t \right) = \dot x\left( t \right) + \dot x\left( t \right) \cdot {u^2}\left( t \right)$ $u\left( t \right) = {\text{arc}}\sin {x\left( t \right)} $ $\dot u\left( t \right) = {{\dot x\left( t \right)} \mathord{\left/ {\vphantom {{\dot x\left( t \right)} {\sqrt {1 - {x^2}\left( t \right)} }}} \right. } {\sqrt {1 - {x^2}\left( t \right)} }}$ $\left. \begin{aligned}& {v\left( t \right) = \sqrt {1 - {x^2}\left( t \right)} } \\ & {w\left( t \right) = {1 / {v\left( t \right)}}} \end{aligned} \right\}$ $\left. \begin{aligned}& {\dot u\left( t \right) = \dot x\left( t \right)w\left( t \right)} \\ & {w\left( t \right)v\left( t \right) - 1 = 0} \\ & {{v^2}\left( t \right) + {x^2}\left( t \right) - 1 = 0} \end{aligned} \right\}$ $u\left( t \right) = {\text{arccos}} {x\left( t \right)} $ $\dot u\left( t \right) = {{ - \dot x\left( t \right)} \mathord{\left/ {\vphantom {{ - \dot x\left( t \right)} {\sqrt {1 - {x^2}\left( t \right)} }}} \right. } {\sqrt {1 - {x^2}\left( t \right)} }}$ $\left.\begin{aligned}& {v\left( t \right) = \sqrt {1 - {x^2}\left( t \right)} } \\ & {w\left( t \right) = - {1 / {v\left( t \right)}}} \end{aligned} \right\}$ $\left. \begin{aligned}& {\dot u\left( t \right) = \dot x\left( t \right)w\left( t \right)} \\ & {w\left( t \right)v\left( t \right) + 1 = 0} \\ & {{v^2}\left( t \right) + {x^2}\left( t \right) - 1 = 0} \end{aligned} \right\}$ $u\left( t \right) = {\text{arctan}}{x\left( t \right)} $ $\dot u\left( t \right) = {{\dot x\left( t \right)} \mathord{\left/ {\vphantom {{\dot x\left( t \right)} {\left( {1 + {x^2}\left( t \right)} \right)}}} \right. } {\left( {1 + {x^2}\left( t \right)} \right)}}$ $\left. \begin{aligned}& {v\left( t \right) = 1 + {x^2}\left( t \right)} \\ & {w\left( t \right) = {1 / {v\left( t \right)}}} \end{aligned}\right\}$ $\left.\begin{aligned}& {\dot u\left( t \right) = \dot x\left( t \right)w\left( t \right)} \\ & {w\left( t \right)v\left( t \right) - 1 = 0} \\ & {v\left( t \right) - {x^2}\left( t \right) - 1 = 0} \end{aligned} \right\}$ -
表 1 RHB法与HDHB法在分岔处的解分布
Table 1 Statistical distribution of solution by the RHB and the HDHB method at bifurcation point
Method Upper
branch/%Lower
branch/%Unstable
branch/%Non-physical/% RHB 57.13 19.71 23.16 0 HDHB 29.14 14.21 15.96 48.69 表 2 各阶RHB-重铸法求解相对论谐振子
Table 2 Solving relativistic harmonic oscillator by the RHB-recast method with different orders
Order of method Amplitude error Computing time/s 25 $4.36 \times {10^{ - 7}}$ 0.73 35 $1.95 \times {10^{ - 9}}$ 0.77 55 $3.87 \times {10^{ - 12}}$ 1.07 表 3 3种方法计算非线性单摆结果对比
Table 3 Comparison of corresponding results of three methods for solving nonlinear pendulum
Methods $M$ Average error Computing time/s RHB-recast 76 1.9 × 10−12 0.70 RHB-Taylor 401 8.9 × 10−10 1.82 HB-AFT 210 1.5 × 10−9 0.69 A1 常见初等超越函数的重铸
A1 Recast form of the most common elementary transcendental functions
Original function Differential relationship Companion variable Recast equation $u\left( t \right) = \exp {x\left( t \right)}$ $\dot u\left( t \right) = \dot x\left( t \right) \cdot u\left( t \right)$ $\dot u\left( t \right) = \dot x\left( t \right) \cdot u\left( t \right)$ $u\left( t \right) = {\log _a} {x\left( t \right)} $ $\dot u\left( t \right) = {{\dot x\left( t \right)} \mathord{\left/ {\vphantom {{\dot x\left( t \right)} {\left( {x\left( t \right) \cdot \ln a} \right)}}} \right. } {\left( {x\left( t \right) \cdot \ln a} \right)}}$ $ v\left( t \right) = {1 \mathord{\left/ {\vphantom {1 {x\left( t \right)}}} \right. } {x\left( t \right)}} $ $\left. \begin{aligned}& {\dot u\left( t \right) = {{\left( {v\left( t \right) \cdot \dot x\left( t \right)} \right)} / {\ln a}}} \\ & {x\left( t \right)v\left( t \right) - 1 = 0} \end{aligned}\right\}$ $\left. \begin{aligned}& {u\left( t \right) = \sin {x\left( t \right)} } \\ & {v\left( t \right) = \cos {x\left( t \right)} } \end{aligned} \right\}$ $\left.\begin{aligned}& {\dot u\left( t \right) = \dot x\left( t \right) \cdot \cos {x\left( t \right)}} \\ & {\dot v\left( t \right) = - \dot x\left( t \right) \cdot \sin {x\left( t \right)} } \end{aligned} \right\}$ $\left. \begin{aligned}& {\dot u\left( t \right) = \dot x\left( t \right) \cdot v\left( t \right)} \\ & {\dot v\left( t \right) = - \dot x\left( t \right) \cdot u\left( t \right)} \end{aligned}\right\}$ $u\left( t \right) = \tan {x\left( t \right)} $ $\dot u\left( t \right) = \dot x\left( t \right) \cdot \left( {1 + {u^2}\left( t \right)} \right)$ $\dot u\left( t \right) = \dot x\left( t \right) + \dot x\left( t \right) \cdot {u^2}\left( t \right)$ $u\left( t \right) = {\text{arc}}\sin {x\left( t \right)} $ $\dot u\left( t \right) = {{\dot x\left( t \right)} \mathord{\left/ {\vphantom {{\dot x\left( t \right)} {\sqrt {1 - {x^2}\left( t \right)} }}} \right. } {\sqrt {1 - {x^2}\left( t \right)} }}$ $\left. \begin{aligned}& {v\left( t \right) = \sqrt {1 - {x^2}\left( t \right)} } \\ & {w\left( t \right) = {1 / {v\left( t \right)}}} \end{aligned} \right\}$ $\left. \begin{aligned}& {\dot u\left( t \right) = \dot x\left( t \right)w\left( t \right)} \\ & {w\left( t \right)v\left( t \right) - 1 = 0} \\ & {{v^2}\left( t \right) + {x^2}\left( t \right) - 1 = 0} \end{aligned} \right\}$ $u\left( t \right) = {\text{arccos}} {x\left( t \right)} $ $\dot u\left( t \right) = {{ - \dot x\left( t \right)} \mathord{\left/ {\vphantom {{ - \dot x\left( t \right)} {\sqrt {1 - {x^2}\left( t \right)} }}} \right. } {\sqrt {1 - {x^2}\left( t \right)} }}$ $\left.\begin{aligned}& {v\left( t \right) = \sqrt {1 - {x^2}\left( t \right)} } \\ & {w\left( t \right) = - {1 / {v\left( t \right)}}} \end{aligned} \right\}$ $\left. \begin{aligned}& {\dot u\left( t \right) = \dot x\left( t \right)w\left( t \right)} \\ & {w\left( t \right)v\left( t \right) + 1 = 0} \\ & {{v^2}\left( t \right) + {x^2}\left( t \right) - 1 = 0} \end{aligned} \right\}$ $u\left( t \right) = {\text{arctan}}{x\left( t \right)} $ $\dot u\left( t \right) = {{\dot x\left( t \right)} \mathord{\left/ {\vphantom {{\dot x\left( t \right)} {\left( {1 + {x^2}\left( t \right)} \right)}}} \right. } {\left( {1 + {x^2}\left( t \right)} \right)}}$ $\left. \begin{aligned}& {v\left( t \right) = 1 + {x^2}\left( t \right)} \\ & {w\left( t \right) = {1 / {v\left( t \right)}}} \end{aligned}\right\}$ $\left.\begin{aligned}& {\dot u\left( t \right) = \dot x\left( t \right)w\left( t \right)} \\ & {w\left( t \right)v\left( t \right) - 1 = 0} \\ & {v\left( t \right) - {x^2}\left( t \right) - 1 = 0} \end{aligned} \right\}$ -
[1] 代洪华, 岳晓奎, 汪雪川. 非线性动力学系统高性能计算方法. 北京: 科学出版社, 2022 (Dai Honghua, Yue Xiaokui, Wang Xuechuan. High-performance Computing Method for Solving Nonlinear Dynamical Systems. Beijing: Science Publishing & Media, 2022 (in Chinese) Dai Honghua, Yue Xiaokui, Wang Xuechuan. High-performance Computing Method for Solving Nonlinear Dynamical Systems. Beijing: Science Publishing & Media, 2022 (in Chinese)
[2] Yan ZP, Dai HH, Wang QS, et al. Harmonic balance methods: A review and recent developments. CMES-Computer Modeling in Engineering & Sciences, 2023, 137(2): 1419-1459
[3] Hall KC, Thomas JP, Clark WS. Computation of unsteady nonlinear flows in cascades using a harmonic balance technique. AIAA Journal, 2002, 40(5): 879-886 doi: 10.2514/2.1754
[4] Liu LP, Thomas JP, Dowell EH, et al. A comparison of classical and high dimensional harmonic balance approaches for a Duffing oscillator. Journal of Computational Physics, 2006, 215(1): 298-320 doi: 10.1016/j.jcp.2005.10.026
[5] LaBryer A, Attar PJ. High dimensional harmonic balance de-aliasing techniques for a Duffing oscillator. Journal of Sound and Vibration, 2009, 324(3-5): 1016-1038 doi: 10.1016/j.jsv.2009.03.005
[6] Dai HH, Yan ZP, Wang XC, et al. Collocation-based harmonic balance framework for highly accurate periodic solution of nonlinear dynamical system. International Journal for Numerical Methods in Engineering, 2023, 124(2): 458-481 doi: 10.1002/nme.7128
[7] Chen Y, Hou L, Chen G, et al. Nonlinear dynamics analysis of a dual-rotor-bearing-casing system based on a modified HB-AFT method. Mechanical Systems and Signal Processing, 2023, 185: 109805
[8] Chen YM, Liu JK, Meng G. An incremental method for limit cycle oscillations of an airfoil with an external store. International Journal of Non-Linear Mechanics, 2012, 47(3): 75-83 doi: 10.1016/j.ijnonlinmec.2011.12.006
[9] Beléndez A, Hernández A, Márquez A, et al. Analytical approximations for the period of a nonlinear pendulum. European Journal of Physics, 2006, 27(3): 539-551 doi: 10.1088/0143-0807/27/3/008
[10] Cameron TM, Griffin JH. An alternating frequency/time domain method for calculating the steady-state response of nonlinear dynamic systems. Journal of Applied Mechanics, 1989, 56(1): 149-154 doi: 10.1115/1.3176036
[11] Krack M, Gross J. Harmonic Balance for Nonlinear Vibration Problems. Cham: Springer International Publishing, 2019: 53-55
[12] Cochelin B, Vergez C. A high order purely frequency-based harmonic balance formulation for continuation of periodic solutions. Journal of Sound and Vibration, 2009, 324(1-2): 243-262 doi: 10.1016/j.jsv.2009.01.054
[13] Karkar S, Cochelin B, Vergez C. A high-order, purely frequency based harmonic balance formulation for continuation of periodic solutions: The case of non-polynomial nonlinearities. Journal of Sound and Vibration, 2013, 332(4): 968-977 doi: 10.1016/j.jsv.2012.09.033
[14] Guillot L, Cochelin B, Vergez C. A generic and efficient Taylor series–based continuation method using a quadratic recast of smooth nonlinear systems. International Journal for Numerical Methods in Engineering, 2019, 119(4): 261-280 doi: 10.1002/nme.6049
[15] Liu G, Lv ZR, Liu JK, et al. Quasi-periodic aeroelastic response analysis of an airfoil with external store by incremental harmonic balance method. International Journal of Non-Linear Mechanics, 2018, 100: 10-19 doi: 10.1016/j.ijnonlinmec.2018.01.004
[16] Huang JL, Zhu WD. A new incremental harmonic balance method with two time scales for quasi-periodic motions of an axially moving beam with internal resonance under single-tone external excitation. Journal of Vibration and Acoustics, 2017, 139(2): 021010
[17] Li H, Ekici K. Supplemental-frequency harmonic balance: A new approach for modeling aperiodic aerodynamic response. Journal of Computational Physics, 2021, 436: 110278
[18] 龚冰清, 郑泽昌, 陈衍茂等. 准周期响应对称破缺分岔点的一种快速计算方法. 力学学报, 2022, 54(11): 3138-3188 (Gong Bingqing, Zheng Zechang, Chen Yanmao, et al. A fast calculation for the symmetry breaking point of quasi-periodic responses. Chinese Journal of Theoretical and Applied Mechanics, 2022, 54(11): 3138-3188 (in Chinese) Gong Bingqing, Zheng Zechang, Chen Yanmao, et al . A fast calculation for the symmetry breaking point of quasi-periodic responses. Chinese Journal of Theoretical and Applied Mechanics,2022 ,54 (11 ):3138 -3188 (in Chinese)[19] Chua L, Ushida A. Algorithms for computing almost periodic steady-state response of nonlinear systems to multiple input frequencies. IEEE Transactions on Circuits and Systems, 1981, 28(10): 953-971 doi: 10.1109/TCS.1981.1084921
[20] Ju R, Fan W, Zhu WD, et al. A modified two-timescale incremental harmonic balance method for steady-state quasi-periodic responses of nonlinear systems. Journal of Computational and Nonlinear Dynamics, 2017, 12(5): 051007
[21] Kim YB, Choi SK. A multiple harmonic balance method for the internal resonant vibration of a non-linear Jeffcott rotor. Journal of Sound and Vibration, 1997, 208(5): 745-761 doi: 10.1006/jsvi.1997.1221
[22] Liu L, Dowell E H, Hall K C. A novel harmonic balance analysis for the Van Der Pol oscillator. International Journal of Non-Linear Mechanics, 2007, 42(1): 2-12 doi: 10.1016/j.ijnonlinmec.2006.09.004
[23] Wang QS, Yan ZP, Dai HH. An efficient multiple harmonic balance method for computing quasi-periodic responses of nonlinear systems. Journal of Sound and Vibration, 2023, 554: 117700
[24] Hall KC, Ekici K, Thomas JP, et al. Harmonic balance methods applied to computational fluid dynamics problems. International Journal of Computational Fluid Dynamics, 2013, 27(2): 52-67 doi: 10.1080/10618562.2012.742512
[25] Dai HH, Schnoor M, Atluri SN. A simple collocation scheme for obtaining the periodic solutions of the duffing equation, and its equivalence to the high dimensional harmonic balance method: subharmonic oscillations. CMES-Computer Modeling in Engineering & Sciences, 2012, 84(5): 459-498
[26] Moreau W, Easther R, Neutze R. Relativistic (an) harmonic oscillator. American Journal of Physics, 1994, 62(6): 531-535 doi: 10.1119/1.17513
[27] Biazar J, Hosami M. An easy trick to a periodic solution of relativistic harmonic oscillator. Journal of the Egyptian Mathematical Society, 2014, 22(1): 45-49 doi: 10.1016/j.joems.2013.04.013
[28] 孟艳平, 孙维鹏, 张皆杰. 相对论谐振子解析逼近解的构造. 吉林大学学报(理学版), 2013, 51(1): 83-88 (Meng Yanping, Sun Weipeng, Zhang Jiejie. Construction of analytical approximate solutions to relativistic harmonic oscillators. Journal of Jilin University (Science Edition), 2013, 51(1): 45-49 (in Chinese) Meng Yanping, Sun Weipeng, Zhang Jiejie. Construction of analytical approximate solutions to relativistic harmonic oscillators. Journal of Jilin University (Science Edition), 2013, 51(1): 45-49 (in Chinese)
[29] Mickens RE. Periodic solutions of the relativistic harmonic oscillator. Journal of Sound and Vibration, 1998, 212(5): 905-908 doi: 10.1006/jsvi.1997.1453
[30] Beléndez A, Pascual C, Méndez DI, et al. Solution of the relativistic (an) harmonic oscillator using the harmonic balance method. Journal of Sound and Vibration, 2008, 311(3-5): 1447-1456 doi: 10.1016/j.jsv.2007.10.010
[31] Zheng ZC, Lv ZR, Chen YM, et al. A modified incremental harmonic balance method combined with Tikhonov regularization for periodic motion of nonlinear system. Journal of Applied Mechanics, 2022, 89(2): 021001
[32] 黄建亮, 张兵许, 陈树辉. 优化迭代步长的两种改进增量谐波平衡法. 力学学报, 2022, 54(5): 1353-1363 (Huang Jianliang, Zhang Bingxu, Chen Shuhui. Two generalized incremental harmonic balance methods with optimization for iteration steps. Chinese Journal of Theoretical and Applied Mechanics, 2022, 54(5): 1353-1363 (in Chinese) Huang Jianliang, Zhang Bingxu, Chen Shuhui . Two generalized incremental harmonic balance methods with optimization for iteration steps. Chinese Journal of Theoretical and Applied Mechanics,2022 ,54 (5 ):1353 -1363 (in Chinese)[33] Thomas J, Dowell EH. Using broyden's method to improve the computational performance of a harmonic balance aeroelastic solution technique//AIAA Scitech 2023 Forum, 2023
[34] Ranganathan A. The Levenberg-Marquardt algorithm. Tutoral on LM Algorithm, 2004, 11(1): 101-110
[35] Liu G, Liu JK, Wang L, et al. Time-domain minimum residual method combined with energy balance for nonlinear conservative systems. Mechanical Systems and Signal Processing, 2022, 170: 108818
[36] Jia QH, Luo Y, Zhou HJ, et al. Nonlinear coupling in an asymmetric pendulum. American Journal of Physics, 2022, 90(8): 594-603 doi: 10.1119/10.0010391
[37] 张哲, 代洪华, 冯浩阳等. 初值约束与两点边值约束轨道动力学方程的快速数值计算方法. 力学学报, 2022, 54(2): 503-516 (Zhang Zhe, Dai Honghua, Feng Haoyang, et al. Efficient numerical method for orbit dynamic functions with initial value and two-point-boundary-value constraints. Chinese Journal of Theoretical and Applied Mechanics, 2022, 54(2): 503-516 (in Chinese) Zhang Zhe, Dai Honghua, Feng Haoyang, et al . Efficient numerical method for orbit dynamic functions with initial value and two-point-boundary-value constraints. Chinese Journal of Theoretical and Applied Mechanics,2022 ,54 (2 ):503 -516 (in Chinese)[38] 王昌涛, 代洪华, 张哲等. 并行加速的局部变分迭代法及其轨道计算应用. 力学学报, 2023, 55(4): 991-1003 (Wang Changtao, Dai Honghua, Zhang Zhe, et al. Parallel accelerated local variational iteration method and its application orbit computation. Chinese Journal of Theoretical and Applied Mechanics, 2023, 55(4): 991-1003 (in Chinese) Wang Changtao, Dai Honghua, Zhang Zhe, et al . Parallel accelerated local variational iteration method and its application orbit computation. Chinese Journal of Theoretical and Applied Mechanics,2023 ,55 (4 ):991 -1003 (in Chinese) -
期刊类型引用(1)
1. 薛皓阳. 硬弹簧Duffing隔振系统的隔振性能与分岔分析. 模具制造. 2024(10): 138-140 . 百度学术
其他类型引用(1)