任意激励下无阻尼系统求解的整体自适应方法_陈东良.pdf

􀪈􀪈􀪈􀪈􀪈􀪈􀪈􀪈􀪈􀪈􀪈􀪈􀪈􀪈􀪈􀪈􀪈􀪈􀪈􀪈􀪈􀪈􀪈􀪈􀪈􀪈􀪈􀪈􀪈􀪈􀪈􀪈􀪈􀪈􀪈􀪈􀪈􀪈􀪈􀪈􀪈􀪈􀪈􀪈􀪈􀪈􀪈􀪈􀪈 􀪈 振 动 与 冲 击 第 ３９ 卷第 ２３ 期ＪＯＵＲＮＡＬ ＯＦ ＶＩＢＲＡＴＩＯＮ ＡＮＤ ＳＨＯＣＫＶｏｌ． ３９ Ｎｏ．２３ ２０２０ 收稿日期 ２０１９ －０５ －０９ 修改稿收到日期 ２０１９ －０７ －１２ 第一作者 陈东良 男，副教授，博士，１９７９ 年生 任意激励下无阻尼系统求解的整体自适应方法 陈东良， 王联甫， 臧 睿， 周长鹤， 宫 臣， 张金东 （哈尔滨工程大学 机电工程学院，哈尔滨 １５０００１） 摘 要针对无阻尼系统求解的计算精度和效率问题，进行了整体自适应精细积分法的研究。 采用 Ｐａｄ 逼近和 Ｓｉｍｐｓｏｎ 公式，分析解方程中矩阵指数和 Ｄｕｈａｍｅｌ 积分项的精细参数之间的关系，据此，提出了振动微分方程组求解的整 体自适应精细积分法。 该方法可以根据不同的精度要求进行自适应积分，对高频激励依旧能保证其计算精度。 这对于求 解高频激励下的振动响应和提高计算效率都是有益的。 数值算例验证了该方法的有效性。 关键词 Ｐａｄ 逼近；Ｄｕｈａｍｅｌ 积分；自适应；精细积分 中图分类号 Ｏ３０２ 文献标志码 ＡＤＯＩ１０． １３４６５／ ｊ． ｃｎｋｉ． ｊｖｓ． ２０２０． ２３． ００７ Ｏｖｅｒａｌｌ ｓｅｌｆ⁃ａｄａｐｔｉｖｅ ｐｒｅｃｉｓｉｏｎ ｉｎｔｅｇｒａｔｉｏｎ ｍｅｔｈｏｄ ｆｏｒ ｓｏｌｖｉｎｇ ａｎ ｕｎｄａｍｐｅｄ ｓｙｓｔｅｍ ｕｎｄｅｒ ａｒｂｉｔｒａｒｙ ｅｘｃｉｔａｔｉｏｎ ＣＨＥＮ Ｄｏｎｇｌｉａｎｇ， ＷＡＮＧ Ｌｉａｎｆｕ， ＺＡＮＧ Ｒｕｉ， ＺＨＯＵ Ｃｈａｎｇｈｅ， ＧＯＮＧ Ｃｈｅｎ， ＺＨＡＮＧ Ｊｉｎｄｏｎｇ （Ｃｏｌｌｅｇｅ ｏｆ Ｍｅｃｈａｎｉｃａｌ ａｎｄ Ｅｌｅｃｔｒｉｃａｌ Ｅｎｇｉｎｅｅｒｉｎｇ， Ｈａｒｂｉｎ Ｅｎｇｉｎｅｅｒｉｎｇ Ｕｎｉｖｅｒｓｉｔｙ， Ｈａｒｂｉｎ １５０００１， Ｃｈｉｎａ） Ａｂｓｔｒａｃｔ Ａｉｍｉｎｇ ａｔ ｃａｌｃｕｌａｔｉｏｎ ｐｒｅｃｉｓｉｏｎ ａｎｄ ｅｆｆｉｃｉｅｎｃｙ ｐｒｏｂｌｅｍｓ ｏｆ ｓｏｌｖｉｎｇ ａｎ ｕｎｄａｍｐｅｄ ｓｙｓｔｅｍ， ｔｈｅ ｏｖｅｒａｌｌ ｓｅｌｆ⁃ ａｄａｐｔｉｖｅ ｐｒｅｃｉｓｉｏｎ ｉｎｔｅｇｒａｔｉｏｎ ｍｅｔｈｏｄ ｗａｓ ｓｔｕｄｉｅｄ． Ｐａｄ ａｐｐｒｏｘｉｍａｔｉｏｎ ａｎｄ Ｓｉｍｐｓｏｎ ｆｏｒｍｕｌａ ｗｅｒｅ ｕｓｅｄ ｔｏ ａｎａｌｙｚｅ ｔｈｅ ｒｅｌａｔｉｏｎｓｈｉｐ ｂｅｔｗｅｅｎ ｍａｔｒｉｘ ｅｘｐｏｎｅｎｔ ａｎｄ ｆｉｎｅ ｐａｒａｍｅｔｅｒｓ ｏｆ Ｄｕｈａｍｅｌ ｉｎｔｅｇｒａｌ ｔｅｒｍ ｉｎ ｓｏｌｖｉｎｇ ｄｉｆｆｅｒｅｎｔｉａｌ ｅｑｕａｔｉｏｎｓ， ａｎｄ ｔｈｅｎ ｔｈｅ ｏｖｅｒａｌｌ ｓｅｌｆ⁃ａｄａｐｔｉｖｅ ｐｒｅｃｉｓｉｏｎ ｉｎｔｅｇｒａｔｉｏｎ ｍｅｔｈｏｄ ｆｏｒ ｓｏｌｖｉｎｇ ａ ｖｉｂｒａｔｉｏｎ ｄｉｆｆｅｒｅｎｔｉａｌ ｅｑｕａｔｉｏｎ ｓｅｔ ｗａｓ ｐｒｏｐｏｓｅｄ． Ｉｔ ｗａｓ ｓｈｏｗｎ ｔｈａｔ ｔｈｉｓｍｅｔｈｏｄ ｃａｎ ｂｅ ｕｓｅｄ ｔｏ ｄｏ ｓｅｌｆ⁃ａｄａｐｔｉｖｅ ｉｎｔｅｇｒａｔｉｏｎ ａｃｃｏｒｄｉｎｇ ｔｏ ｄｉｆｆｅｒｅｎｔ ｐｒｅｃｉｓｉｏｎ ｒｅｑｕｉｒｅｍｅｎｔｓ， ａｎｄ ｋｅｅｐ ｉｔｓ ｃｏｍｐｕｔａｔｉｏｎａｌ ａｃｃｕｒａｃｙ ｕｎｄｅｒ ｈｉｇｈ ｆｒｅｑｕｅｎｃｙ ｅｘｃｉｔａｔｉｏｎ； ｔｈｉｓ ｍｅｔｈｏｄ ｉｓ ｂｅｎｅｆｉｃｉａｌ ｔｏ ｓｏｌｖｅ ｖｉｂｒａｔｉｏｎ ｒｅｓｐｏｎｓｅｓ ｕｎｄｅｒ ｈｉｇｈ ｆｒｅｑｕｅｎｃｙ ｅｘｃｉｔａｔｉｏｎ ａｎｄ ｉｍｐｒｏｖｅ ｃｏｍｐｕｔａｔｉｏｎａｌ ｅｆｆｉｃｉｅｎｃｙ． Ｎｕｍｅｒｉｃａｌ ｅｘａｍｐｌｅｓ ｖｅｒｉｆｉｅｄ ｔｈｅ ｅｆｆｅｃｔｉｖｅｎｅｓｓ ｏｆ ｔｈｅ ｐｒｏｐｏｓｅｄ ｍｅｔｈｏｄ． Ｋｅｙ ｗｏｒｄｓ Ｐａｄ ａｐｐｒｏｘｉｍａｔｉｏｎ； Ｄｕｈａｍｅｌ ｉｎｔｅｇｒａｌ； ｓｅｌｆ⁃ａｄａｐｔｉｖｅ； ｐｒｅｃｉｓｉｏｎ ｉｎｔｅｇｒａｔｉｏｎ 结构动力响应、受激振动以及优化控制等工程领 域的许多瞬态历程问题，都可以直接，或经过一定变换 转化为一阶线性常系数微分方程组。 而矩阵指数及 Ｄｕｈａｍｅｌ 积分项的精细计算是其求解的关键，在计算过 程中，精度要求与激励的频率都会对精细参数的选择 产生影响。 钟万勰［１］提出的精细积分法解决了矩阵指数的精 细计算问题，并得到了广泛的应用。 许多文献研究了 使用 Ｐａｄ 逼近对矩阵指数展开的优越性，并进行了误 差分析［２⁃４］。 文献［５］基于误差分析公式提出了矩阵指 数计算时的参数自适应选择方法，具有良好的效果。 对于 Ｄｕｈａｍｅｌ 积分项，一些文献采用不同的插值方法 进行运算［６⁃７］，但这些方法无法适应高频激励，或需划 分很小的积分区间并以很大的计算量为代价。 文献 ［８］使用 Ｒｏｍｂｅｒｇ 积分法进行积分项的自适应求解，该 方法需事先指定一个确定的积分步长，其仅能控制单 个积分区间的精度范围，由于传递误差的存在，后续节 点解的精度无法得到有效保证。 本文采用基于 Ｓｉｍｐｓｏｎ 公式的自适应积分法，对 Ｄｕｈａｍｅｌ 积分项进行自适应计算。 在自适应积分得到 的积分子区间上，利用谭述君等的研究的方法确定矩 阵指数计算的精细参数，并基于 Ｐａｄ 逼近对矩阵指数 进行精细计算，建立各精细参数之间的关系，从而确定 整体的自适应积分方法，降低计算量，提高计算效率， 实现了给定容差下的自适应精细计算。 １ 整体自适应精细积分法 无阻尼系统方程、振动方程等瞬态历程问题都可 以通过引入状态向量的方式转化为一阶线性常系数微 分方程组 ｖ （ｔ） ＝ Ａｖ（ｔ） ＋ ｒ（ｔ） ｖ（０） ＝ ｖ０ { （１） 式中ｒ（ｔ）是引入的状态向量，Ａ 是系数矩阵，ｒ（ｔ）是外 力载荷或激励向量，ｖ０是状态向量的初值。 其解为包 含 Ｄｕｈａｍｅｌ 积分项的齐次解与特解的组合 ｖ（ｔ） ＝ ｅＡ（ｔ－ｔ０）ｖ０＋ ｅＡｔ∫ ｔ ｔ０ｅ －Ａτｆ（τ）ｄτ （２） １． １ 基于 Ｐａｄ 逼近的矩阵指数精细计算 Ｐａｄ 级数逼近的整体误差要小于多项式逼近，对 矩阵指数的（ｐ，ｑ）阶 Ｐａｄ 逼近定义为［１］ Ｒｐｑ（Ａ） ＝ ［Ｄｐｑ（Ａ）］ －１Ｎ ｐｑ（Ａ） （３） 式中， Ｎｐｑ（Ａ） ＝∑ ｐ ｊ ＝０ （ｐ ＋ ｑ － ｊ）ｐ （ｐ ＋ ｑ）ｊ（ｐ － ｊ）Ａ ｊ （４） Ｄｐｑ（Ａ） ＝∑ ｑ ｊ ＝０ （ｐ ＋ ｑ － ｊ）ｑ （ｐ ＋ ｑ）ｊ（ｑ － ｊ）（ － Ａ） ｊ （５） Ｐａｄ 逼近在原点附近的误差最小，因此，为保证逼 近效果，使用加权平方法来保证这一点 Ｆｐｑ＝Ｒｐｑ Ａ ２Ｎ ２Ｎ ＝ （Ｒｐｑ（Ａ′））２Ｎ（６） 式中Ａ′ ＝ Ａ ２Ｎ，Ｎ 为加权数，也称精细参数。 当 ｐ ＝ ｑ 时，可以取得更好的精度效果，因此本文令 ｐ ＝ ｑ。 计算 时需将 Ｆｐｑ的主部与增量分离，避免在作乘方运算时精 度丧失。 然后使用钟万勰院士提出的方法精细计算。 Ｍｏｌｅｒ 等对使用 Ｐａｄ 逼近计算矩阵指数时的误差 给出了如下证明 当 Ａ ∞／２ Ｎ≤１／２ 时，存在矩阵 Ｅ，有 Ｅ ∞ ≤ ε（Ｎ，ｑ）Ａ ∞ （７） 式中ε（Ｎ，ｑ） ＝８ Ａ ∞ ２Ｎ ２ｑ （ｑ）２ （２ｑ） （２ｑ ＋１） （８） 谭述君等给出了精细参数（Ｎ，ｑ）自适应选择的递 推过程 实际的计算中容差一般选取为较小的数，因此误 差界取 ε（Ｎ，ｑ） Ａ ∞即可满足精度要求。 若给定容差 为 ｔｏｌ，则 Ｅｒｒ（Ｎ，ｑ） ＝ ε（Ｎ，ｑ） ＝ ε（Ｎ，ｑ）Ａ ∞ ≤ ｔｏｌ （９） 然后结合以下三组递推关系式即可实现（Ｎ，ｑ）的 自动选择 Ｅｒｒ（Ｎ＋１，ｑ） Ｅｒｒ（Ｎ，ｑ） ＝ε（Ｎ＋１，ｑ） ε（Ｎ，ｑ） ＝ １ ２２ｑ ＝ρＮ（Ｎ，ｑ） Ｅｒｒ（Ｎ，ｑ ＋ １） Ｅｒｒ（Ｎ，ｑ） ＝ ε（Ｎ，ｑ ＋ １） ε（Ｎ，ｑ） ＝ Ａ ∞ ２Ｎ ２ １ ４（２ｑ ＋ １）（２ｑ ＋ ３） ＝ ρｑ（Ｎ，ｑ） （１０） ρＮ（Ｎ，ｑ ＋ １） ＝ １ ４ ρＮ（Ｎ，ｑ） ρｑ（Ｎ，ｑ ＋ １） ＝ ２ｑ ＋ １ ２ｑ ＋ ５ρｑ（Ｎ，ｑ） （１１） ρＮ（Ｎ ＋ １，ｑ） ＝ ρＮ（Ｎ，ｑ） ρｑ（Ｎ ＋ １，ｑ） ＝ １ ４ ρｑ（Ｎ，ｑ） { （１２） 初始化参数 ｑ０＝１，Ｎ０＝ ｍａｘ（０，ｌｏｇ２Ａ ∞ ＋１）。 在 编程计算时，需将 Ｎ０向上圆整为一个整数。 １． ２ 基于 Ｓｉｍｐｓｏｎ 公式的自适应积分 当激励曲线变化剧烈时，需使用更小的积分步长 来达到给定的精度要求，而在变化相对平缓的区段，则 可使用较大的积分步长，这相对使用等距节点的积分 更有效率。 据此，文献［９］使用辛普森公式提出了一种 自适应积分的方法 辛普森公式使用［ａｋ，ｂｋ］上的两个子区间表述为 Ｓ（ａｋ，ｂｋ） ＝ ｈ ３ （ｆ（ａｋ） ＋ ４ｆ（ｃｋ） ＋ ｆ（ｂｋ））（１３） 式中ｃｋ是区间［ａｋ，ｂｋ］的中心，ｈ ＝ （ｂｋ－ ａｋ） ／２。 ［ａｋ，ｂｋ］上 ４ 个子区间的组合 Ｓｉｍｐｓｏｎ 公式可以先 将区间划分为两个相等的子区间［ａｋ１，ｂｋ１］ 和［ａｋ２， ｂｋ２］，然后在两段上递归应用式（１３）。 Ｓ（ａｋ１，ｂｋ１） ＋ Ｓ（ａｋ２，ｂｋ２） ＝ ｈ ６ （ｆ（ａｋ１） ＋ ４ｆ（ｃｋ１） ＋ ｆ（ｂｋ１）） ＋ ｈ ６ （ｆ（ａｋ２） ＋ ４ｆ（ｃｋ２） ＋ ｆ（ｂｋ２））（１４） 式中ａｋ１＝ ａｋ，ｂｋ１＝ ａｋ２＝ ｃｋ，ｂｋ２＝ ｂｋ，ｃｋ１是［ａｋ１，ｂｋ１］的中 点，ｃｋ２是［ａｋ２，ｂｋ２］的中点。 对于精度测试，指定［ａｋ，ｂｋ］上的容差，若有 １ １０ Ｓ（ａｋ１，ｂｋ１） ＋ Ｓ（ａｋ２，ｂｋ２） － Ｓ（ａｋ，ｂｋ） ＜ ｔｏｌ（１５） 则认为 ∫ ｂｋ ａｋｆ（ｘ）ｄｘ － Ｓ（ａｋ１，ｂｋ１） － Ｓ（ａｋ２，ｂｋ２） ＜ ｔｏｌ（１６） 通过使用式（１３）和式（１４）实现自适应积分。 从 ｛［ａ０，ｂ０］，ｔｏｌ｝开始，将该区间划分为两个子区间，分别 为［ａ０１，ｂ０１］和［ａ０２，ｂ０２］。 如果能够通过精度测试式 （１５），则在区间［ａ０，ｂ０］上使用式（１４），过程结束。 如 果未通过测试，则将［ａ０１，ｂ０１］和［ａ０２，ｂ０２］分别标记为 ［ａ１，ｂ１］和［ａ２，ｂ２］，对应的容差分别变为 ｔｏｌ１＝ ｔｏｌ２＝ １ ２ ｔｏｌ，然后重复上述过程。 １． ３ 整体自适应求解方法 对解方程（２）应用整体自适应精细积分时，需确定 矩阵指数和 Ｄｕｈａｍｅｌ 积分之间的精细参数对应关系， 从而保证整体的容差要求。 其过程可分为以下三步 （１） 基于 Ｓｉｍｐｓｏｎ 公式计算 Ｄｕｈａｍｅｌ 项自适应积 分，得到划分的子区间及对应的区间积分和容差。 其 中 Ｄｕｈａｍｅｌ 积分项中的矩阵指数使用当前子区间的容 差确定精细参数（Ｎ，ｑ）。 ８４振 动 与 冲 击 ２０２０ 年第 ３９ 卷 （２） 根据第一步计算结果，从第一个子区间开始， 计算 ｅＡ（ｔ － ｔ０）与 ｅＡｔ，其中，（ｔ － ｔ０）和 ｔ 分别为当前子区间 的步长和右端点。 ［Ｎ，ｑ］计算时容差的选取与当前子 区间保持一致。 （３） 将前两步中所得结果代入式（２）得到当前子 区间的节点解，然后以当前节点解为初值计算下一个 子区间，重复此过程，直至获取所有子区间上的节点 解。 计算流程，见图 １。 图 １ 整体自适应积分计算流程图 Ｆｉｇ． １ Ｆｌｏｗ ｃｈａｒｔ ｏｆ ｏｖｅｒａｌｌ ｓｅｌｆ⁃ａｄａｐｔｉｖｅ ｉｎｔｅｇｒａｌ ｍｅｔｈｏｄ ２ 数值算例 ２． １ 算例 １ 考虑一个外激励为正余弦函数组合的简单系统， 其初值问题 ２０ ０１ [] ｘ １ ｘ ２ { }＋ ６－ ２ － ２４ [] ｘ１ ｘ２ { } ＝ ０ ２ [ ]ｓｉｎ ωｔ ＋ ２ ０ [ ]ｃｏｓ ωｔ ｘ１（０） ＝ ｘ２（０） ＝ ０ ｘ １（０） ＝ ｘ ２（０） ＝ ０ ω ＝ ２πｆ （１７） 式中ｆ 为外激励的频率。 为简化计算，引入状态向量 ｖ ＝ ｛ｘ１（ｔ） ｘ２（ｔ） ｘ １（ｔ） ｘ ２（ｔ）｝后，式（１７）化为一 阶常系数非其次微分方程组 ｄ ｄｔ ｘ１（ｔ） ｘ２（ｔ） ｘ １（ｔ） ｘ ２（ｔ） ＝ ００１０ ０００１ － ３１００ ２－ ４００ ｘ１（ｔ） ｘ２（ｔ） ｘ １（ｔ） ｘ ２（ｔ） ＋ ０ ０ ０ ２ ｓｉｎ ωｔ ＋ ０ ０ １ ０ ｃｏｓ ωｔ（１８） 为验证本文方法计算的准确性，需找出精确解为 参照进行对比。 当外激励为正／ 余弦函数组合时有解 析解［１０］ ｖｋ＋１＝ Ｔｖｋ－ Ｔ［ａｓｉｎ（ωｔｋ） ＋ ｂｃｏｓ（ωｔｋ）］ ＋ ａｓｉｎ（ωｔｋ＋１） ＋ ｂｃｏｓ（ωｔｋ＋１） ａ ＝ （ω２Ｉ ＋ Ａ２） －１（ωｒ ２ － Ａｒ１） ｂ ＝ （ω２Ｉ ＋ Ａ２） －１（ － ωｒ １ － Ａｒ２） （１９） 式中Ｔ ＝ ｅＡ（ｔ － ｔ０），ｒ１和 ｒ２分别是正弦和余弦项的系数 矩阵。 本例以式（１９）作为精确解参照。 给定积分区间为［０，１］，容差 ｔｏｌ ＝０． ０００ ０００ ００１， 分别计算 ｆ ＝ ［１ ５ １０ １５］时的积分结果。 取精确 解和本文解在 ｔ ＝０． ５ 和 ｔ ＝１ 的两个值进行比较，结果 如表 １ 所示。 本例主要验证该方法对于激励频率的自适应性。 在实际应用中，外部激励存在任意性，且可能出现瞬时 高频或持续高频现象。 因此本例采用变化频率可控的 简谐函数组合作为激励源，通过控制频率参数 ｆ 达到验 证目的。 以上结果表明，对于不同频率的激振函数，该 方法均可根据给定的时程区间和容差有效划分子区 间，从而达到精度要求。 从两个节点解的结果可知，对于不同频率的激励 曲线，其均可达到小数点后 ９ 位及以上精度。 虽未与 给定的 ９ 位容差要求完全一致，但这并非表明精度设 计失败了。 由于给定的容差是整个积分区间的总体误 差之和，因此随着区间的细分，单个子区间上的误差必 然小于给定的容差。 从表 １ 的结果还可看出，随着频率的增加，有效精 度的位数也在增加。 这是由于频率越高，激励曲线变 化越激烈，自适应计算需划分的子区间也就越多。 因 此，给定的容差在这种情况下被分割为更小的区间容 差，所以单个子区间精度随频率增大而有所提高。 ９４第 ２３ 期陈东良等 任意激励下无阻尼系统求解的整体自适应方法 表 １ 不同激励频率下的节点解 Ｔａｂ． １ Ｎｏｄｅ ｓｏｌｕｔｉｏｎ ａｔ ｄｉｆｆｅｒｅｎｔ ｅｘｃｉｔａｔｉｏｎ ｆｒｅｑｕｅｎｃｉｅｓ ｔ ＝０． ５ｔ ＝１ Ｅｘａｃｔ ｓｏｌｕｔｉｏｎＴｈｉｓ ＰａｐｅｒＥｘａｃｔ ｓｏｌｕｔｉｏｎＴｈｉｓ ｐａｐｅｒ ｆ ＝１ ０． ０４７ ６９１ ５３１ ０６８ ２００． ０４７ ６９１ ５３１ ０５１ ６５０． ００４ ６６７ ８０７ １４６ １５０． ００４ ６６７ ８０７ １６１ ２８ ０． １５２ ５７３ ２６９ ３５８ ８４０． １５２ ５７３ ２６９ ３７８ ８９０． １７９ ８１１ ７６１ １６６ ６８０． １７９ ８１１ ７６１ １８５ １４ ｆ ＝５ ０． ００２ ８６９ ７２８ １４３ ５００． ００２ ８６９ ７２８ １４２ ３００． ００ ６２２ ７１０ ３２５ ９９２０． ００６ ２２７ １０３ ２６３ ２１ ０． ０２７ １３８ １０５ ６９６ ０２０． ０２７ １３８ １０５ ７０４ １７０． ０３０ ３８７ ２８２ ６９６ ０００． ０３０ ３８７ ２８２ ７０４ ７７ ｆ ＝１０ ５． １６０ ０３３ ３２３ ３ １０ －４ ５． １６０ ０３３ ３２６ ０ １０ －４ ０． ００３ ３９７ ３４６ ６５４ ２５０． ００３ ３９７ ３４６ ６５５ ６３ ０． ０１３ ４７５ ８４９ ９９１ １３０． ０１３ ４７５ ８４９ ９９５ ２１０． ０１５ ０２１ ２９９ ８２６ ２３０． ０１５ ０２１ ２９９ ８３０ ７０ ｆ ＝１５ ５． ８３８ ５６８ ０２３ ０ １０ －４ ５． ８３８ ５６８ ０２２ ８ １０ －４ ０． ００２ ３２７ ３９５ ３２１ ３９０． ００２ ３２７ ３９５ ３２１ ７６ ０． ００８ ９０２ ７０６ ０４４ ６５０． ００８ ９０２ ７０６ ０４５ ５２０． ００９ ９８０ １５４ ９７７ ２４０． ００９ ９８０ １５４ ９７８ ２２ 取 ｆ ＝１５ 时的 ２ ０４８ 个节点误差结果绘制 ｘ１和 ｘ２ 的误差曲线，如图 ２ 所示。 图 ２ 相对误差曲线 Ｆｉｇ． ２ Ｒｅｌａｔｉｖｅ ｅｒｒｏｒ ｃｕｒｖｅ 从误差曲线中可以看出，随着节点数的增加，误差 曲线呈一定周期性的波动。 这是由于非齐次项为周期 函数，当激励曲线变化剧烈时，误差也相应地增大，因 此也呈周期性变化，但其误差仍然控制在给定的容差 范围之内。 从原理上讲，误差会随着节点的传递不断 增大，但从图中可知，误差传递对其的影响远小于给定 的容差范围，因此认定该方法是有效的。 ２． ２ 算例 ２ 如图 ３ 所示的三自由度无阻尼弹簧质量系统，ｘ１， ｘ２，ｘ３表示该系统稳态响应。 由于本文主要讨论系统 图 ３ 三自由度弹簧质量系统 Ｆｉｇ． ３ Ｔｈｒｅｅ⁃ｄｅｇｒｅｅｓ⁃ｏｆ⁃ｆｒｅｅｄｏｍ ｓｐｒｉｎｇ⁃ｍａｓｓ ｓｙｓｔｅｍ 求解方法的自适应性，因此在这里不考虑其达到稳态 之前的过程。 该无阻尼系统的强迫振动方程如式（２０） 所示 Ｍｘ ＋ Ｋｘ ＝ Ｐ（ｔ）（２０） 式中，Ｐ（ｔ） ＝ ［Ｐ１（ｔ） Ｐ２（ｔ） Ｐ３（ｔ）］为激振矢量。 取 ｍ ＝２，ｋ ＝１，Ｐ１（ｔ） ＝ １３（ｔ － ｔ２）ｅ －３ｔ／２，Ｐ ２（ｔ） ＝ Ｐ３（ｔ） ＝ ０，则式（２０）可以表示为 ２００ ０２０ ００２ ｘ １ ｘ ２ ｘ ３ ＋ ３－ １０ － １２－ １ ０－ １３ ｘ１ ｘ２ ｘ３ ＝ １３（ｔ － ｔ２）ｅ －３ｔ／２ ０ ０ （２１） 用 Ｐ１（ｔ） ＝１３（ｔ － ｔ２）ｅ －３ｔ／２来模拟该系统受到的任 意激励，该函数在 ｔ ＝ ［０，４］的曲线，如图 ４ 所示。 图 ４ 任意激励曲线 Ｆｉｇ． ４ Ａｒｂｉｔｒａｒｙ ｅｘｃｉｔａｔｉｏｎ ｃｕｒｖｅ 引入状态向量 ｖ ＝ ｛ｘ１ ｘ２ ｘ３ ｘ １ ｘ ２ ｘ ３｝，将 式（２１）化为 ｄ ｄｔ ｘ１ ｘ２ ｘ３ ｘ １ ｘ ２ ｘ ３ ＝ ０００１００ ００００１０ ０００００１ － ３ ２ １ ２ ００００ １ ２ － １ １ ２ ０００ ０ １ ２ － ３ ２ ０００ ｘ１ ｘ２ ｘ３ ｘ １ ｘ ２ ｘ ３ ＋ ０ ０ ０ １３ ２ （ｔ － ｔ２）ｅ －３ｔ／２ ０ ０ （２２） ０５振 动 与 冲 击 ２０２０ 年第 ３９ 卷 式（２２）即 ｖ ＝ Ａｖ ＋ ｒ（ｔ），其中，Ａ 为 ６ ６ 系数矩阵， ｒ（ｔ） ＝ ０ ０ ０ １３ ２ （ｔ － ｔ２）ｅ － ３ｔ ２ ０ ０ [] Ｔ 。 除本文方法外，矩阵指数还可使用特征值法求解 假设 Ａ 为三阶矩阵，［λ１ λ２ λ３］和 Ｐ ＝ ｛η１，η２， η３｝分别为其一组对应的特征值与特征向量，则有 ｅＡ＝ Ｐｄｉａｇ（ｅλ１，ｅλ２，ｅλ３）Ｐ －１ （２３） 使用特征值解法求出式（２）中 ｅ － Ａτ，然后求积分 Ｓ（ｔ） ＝∫ｅ －Ａｔｒ（ｔ）ｄｔ （２４） 式中，Ｓ（ｔ） ＝ ［Ｓ１（ｔ） Ｓ２（ｔ） Ｓ３（ｔ） Ｓ４（ｔ） Ｓ５（ｔ） Ｓ６（ｔ）］。 本例取初值 ｖ０＝ ｛０｝，则该算例在时程［ａ，ｂ］ 上的精确解为 ｖａｂ＝ ｅＡｂ（Ｓ（ｂ） － Ｓ（ａ））（２５） 为保证式（２５）的结果为精确解，ｅＡｂ仍选用特征值 解法。 现给定 ｔ ＝ ［０，４］，容差 ｔｏｌ ＝ ０． ０００ ００１，分别使 用本文的自适应计算方法和精确解式（２５）计算，所得 结果划分子区间为 ３８ 个，对应 Ｓ（ｔ）的划分结果如图 ５ 所示。 （ａ） （ｂ） （ｃ） （ｄ） （ｅ） （ｆ） 图 ５ 子区间划分结果 Ｆｉｇ． ５ Ｓｕｂｉｎｔｅｒｖａｌ ｄｉｖｉｓｉｏｎ ｒｅｓｕｌｔ 从子区间划分结果可以看出，有些较平缓的曲线 段依旧会划分较多的子区间。 这是由于本方法保证的 是在相同子区间下 ６ 条 Ｓ（ｔ）曲线的最大误差满足给定 的容差要求，因此只要该子区间任意一条曲线未达到 精度要求，就需要作区间细分。 这同时也保证 Ｓ（ｔ）的 ６ 条曲线使用共同的划分结果，方便数据处理。 本例主 要为了说明整体自适应方法对任意激励的敏感性，因 此不再具体列出各节点的计算结果。 图 ５ 表明该方法 对激励曲线具有良好的适应性。 将自适应计算的结果与精确解对比得误差曲线， 如图 ６ 所示。 图 ６ 绝对误差曲线 Ｆｉｇ． ６ Ａｂｓｏｌｕｔｅ ｅｒｒｏｒ ｃｕｒｖｅ 从图 ６ 的误差曲线可知，三个未知量的单节点误 差依然远小于给定的容差，满足精度要求。 ３ 结 论 本文介绍了矩阵指数的自适应精细计算和 Ｄｕｈａｍｅｌ 项的自适应积分，对于任意激励下无阻尼系统响应的初 值问题，讨论了其解方程中矩阵指数与积分项之间自适 应参数之间的关系，并给出了整体自适应积分的步骤。 数值算例计算了在不同激励频率下该方法的自适应性精 度及误差变化，并给出了在特定初始条件下子区间划分 结果。 算例表明该方法能够对任意激励进行自适应区间 划分，满足给定的容差要求，并可应用于较高的频率范 围。 由于计算过程中的精细参数均为自适应选择，因此 该方法同时也提升了计算的效率。 参 考 文 献 ［ １］ 钟万勰． 结构动力方程的精细时程积分法［Ｊ］． 大连理工 大学学报， １９９４，３４（２）１３１⁃１３６． ＺＨＯＮＧ Ｗａｎｘｉｅ．Ｏｎ ｐｒｅｃｉｓｅ ｔｉｍｅ⁃ｉｎｔｅｇｒａｔｉｏｎ ｍｅｔｈｏｄ ｆｏｒ ｓｔｒｕｃｔｕｒａｌ ｄｙｎａｍｉｃｓ ［ Ｊ］．Ｊｏｕｒｎａｌ ｏｆ Ｄａｌｉａｎ Ｕｎｉｖｅｒｓｉｔｙ ｏｆ Ｔｅｃｈｎｏｌｏｇｙ， １９９４， ３４（２） １３１⁃１３６． ［ ２］ 刘勇， 沈为平． 精细时程积分中状态转换矩阵的自适应算 法［Ｊ］． 振动与冲击， １９９５，１４（２）８２⁃８５． ＬＩＵ Ｙｏｎｇ， ＳＨＥＮ Ｗｅｉｐｉｎｇ．Ａｄａｐｔｉｖｅ ａｌｇｏｒｉｔｈｍ ｆｏｒ ｓｔａｔｅ⁃ ｔｒａｎｓｉｔｉｏｎ ｍａｔｒｉｘ ｉｎ ｐｒｅｃｉｓｅ ｔｉｍｅ ｓｔｅｐ ｉｎｔｅｇｒａｔｉｏｎ ｍｅｔｈｏｄ［Ｊ］． Ｊｏｕｒｎａｌ ｏｆ Ｖｉｂｒａｔｉｏｎ ａｎｄ Ｓｈｏｃｋ， １９９５，１４（２） ８２⁃８５． （下转第 ７０ 页） １５第 ２３ 期陈东良等 任意激励下无阻尼系统求解的整体自适应方法 ｒｅｓｏｎａｎｃｅ ｉｎ ｔｈｅ ｇｒｏｗｔｈ ｏｆ ａ ｔｕｍｏｒ ｉｎｄｕｃｅｄ ｂｙ ｃｏｒｒｅｌａｔｅｄ ｎｏｉｓｅｓ ［ Ｊ］．Ｃｈｉｎｅｓｅ Ｓｃｉｅｎｃｅ Ｂｕｌｌｅｔｉｎ， ２００５， ５０ （ ２０ ） ２４７５⁃２４７８． ［１１］ ＭＡＳＯＮ Ｊ， ＬＩＮＤＮＥＲ Ｊ Ｆ， ＮＥＦＦ Ｊ， ｅｔ ａｌ． Ｐｕｌｓｅ⁃ｅｎｈａｎｃｅｄ ｓｔｏｃｈａｓｔｉｃ ｒｅｓｏｎａｎｃｅ［Ｊ］． Ｐｈｙｓｉｃｓ Ｌｅｔｔｅｒｓ Ａ， ２０００， ２７７（１） １３⁃１７． ［１２］ ＧＵＯ Ｆ，ＬＵＯ Ｘ Ｄ，ＬＩ Ｓ Ｆ，ｅｔ ａｌ． Ｓｔｏｃｈａｓｔｉｃ ｒｅｓｏｎａｎｃｅ ｉｎ ａ ｍｏｎｏｓｔａｂｌｅｓｙｓｔｅｍｄｒｉｖｅｎｂｙｓｑｕａｒｅ⁃ｗａｖｅｓｉｇｎａｌａｎｄ ｄｉｃｈｏｔｏｍｏｕｓ ｎｏｉｓｅ［Ｊ］． Ｃｈｉｎｅｓｅ Ｐｈｙｓｉｃｓ Ｂ， ２０１０， １９（８） １５６⁃１６２． ［１３］ ＸＵ Ｗ，ＪＩＮ Ｙ Ｆ，ＬＩ Ｗ，ｅｔ ａｌ． Ｓｔｏｃｈａｓｔｉｃ ｒｅｓｏｎａｎｃｅ ｉｎ ａｎ ａｓｙｍｍｅｔｒｉｃ ｂｉｓｔａｂｌｅ ｓｙｓｔｅｍ ｄｒｉｖｅｎ ｂｙ ｍｕｌｔｉｐｌｉｃａｔｉｖｅ ａｎｄ ａｄｄｉｔｉｖｅ ｎｏｉｓｅ ［ Ｊ ］．Ｃｈｉｎｅｓｅ Ｐｈｙｓｉｃｓ， ２００５， １４ （ ６ ） １０７７⁃１０８１． ［１４］ ＪＩＮ Ｘ Ｆ， ＸＵ Ｗ， ＸＵ Ｍ．Ｓｔｏｃｈａｓｔｉｃ ｒｅｓｏｎａｎｃｅ ｉｎ ａｎ ａｓｙｍｍｅｔｒｉｃ ｂｉｓｔａｂｌｅｓｙｓｔｅｍｓｕｂｊｅｃｔｔｏｆｒｅｑｕｅｎｃｙｍｉｘｉｎｇ ｐｅｒｉｏｄｉｃ ｆｏｒｃｅ ａｎｄ ｎｏｉｓｅ［Ｊ］． Ｃｈｉｎｅｓｅ Ｐｈｙｓｉｃｓ Ｌｅｔｔｅｒｓ， ２００５， ２２（５）１０６１⁃１０６４． ［１５］ 董小娟． 含关联噪声与时滞项的非对称双稳系统的随机 共振［Ｊ］． 物理学报， ２００７， ５６（１０）５６１８⁃５６２２． ＤＯＮＧ Ｘｉａｏｊｕａｎ．Ｓｔｏｃｈａｓｔｉｃ ｒｅｓｏｎａｎｃｅ ｉｎ ａｎ ａｓｙｍｍｅｔｒｉｃ ｂｉｓｔａｂｌｅ ｓｙｓｔｅｍ ｗｉｔｈ ｔｉｍｅ⁃ｄｅｌａｙｅｄ ｆｅｅｄｂａｃｋ ａｎｄ ｃｏｒｒｅｌａｔｅｄ ｎｏｉｓｅｓ［Ｊ］． Ａｃｔａ Ｐｈｙｓｉｃａ Ｓｉｎｉｃａ， ２００７， ５６（１０）５６１８⁃５６２２． ［１６］ ＺＨＡＮＧ Ｌｅｉ， ＺＨＥＮＧ Ｗｅｎｂｉｎ， ＦＥＩ Ｘｉｅ， ｅｔ ａｌ． Ｅｆｆｅｃｔ ｏｆ ｔｈｅ ｃｏｒｒｅｌａｔｉｏｎ ｂｅｔｗｅｅｎ ｉｎｔｅｒｎａｌ ｎｏｉｓｅ ａｎｄ ｅｘｔｅｒｎａｌ ｎｏｉｓｅ ｏｎ ｌｏｇｉｃａｌ ｓｔｏｃｈａｓｔｉｃ ｒｅｓｏｎａｎｃｅ ｉｎ ｂｉｓｔａｂｌｅ ｓｙｓｔｅｍｓ［Ｊ］． Ｐｈｙｓｉｃａｌ Ｒｅｖｉｅｗ Ｅ，２０１７，９６（５） ０５２２０３． ［１７］ ＸＩＥ Ｃ Ｗ， ＭＥＩ Ｄ Ｃ， ＣＡＯ Ｌ， ｅｔ ａｌ． Ｅｓｃａｐｅ ｏｆ ｂｒｏｗｎｉａｎ ｐａｒｔｉｃｌｅｓ ｉｎ ａｎ ａｓｙｍｍｅｔｒｉｃ ｂｉｓｔａｂｌｅ ｓａｗｔｏｏｔｈ ｐｏｔｅｎｔｉａｌ ｄｒｉｖｅｎ ｂｙ ｃｒｏｓｓ⁃ｃｏｒｒｅｌａｔｅｄ ｎｏｉｓｅｓ［Ｊ］． Ｅｕｒｏｐｅａｎ Ｐｈｙｓｉｃａｌ Ｊｏｕｒｎａｌ Ｂ， ２００３， ３３（１）８３⁃８６． ［１８］ 周丙常，徐伟． 周期混合信号和噪声联合激励下的非对称 双稳 系 统 的 随 机 共 振 ［ Ｊ］． 物 理 学 报， ２００７ （ １０ ） ５６２３⁃５６２８． ＺＨＯＵ Ｂｉｎｇｃｈａｎｇ， ＸＵ Ｗｅｉ．Ｓｔｏｃｈａｓｔｉｃ ｒｅｓｏｎａｎｃｅ ｉｎ ａｎ ａｓｙｍｍｅｔｒｉｃ ｂｉｓｔａｂｌｅ ｓｙｓｔｅｍ ｄｒｉｖｅｎ ｂｙ ｍｉｘｅｄ ｐｅｒｉｏｄｉｃ ｆｏｒｃｅ ａｎｄ ｎｏｉｓｅｓ［Ｊ］． Ａｃｔａ Ｐｈｙｓｉｃａ Ｓｉｎｉｃａ， ２００７（１０）５６２３⁃５６２８． ［１９］ 焦尚彬， 任超， 黄伟超， 等． α 稳定噪声环境下多频微弱 信号检测的参数诱导随机共振现象［Ｊ］． 物理学报， ２０１３， ６２（２１） ３４⁃４３． ＪＩＡＯ Ｓｈａｎｇｂｉｎ， ＲＥＮ Ｃｈａｏ， ＨＵＡＮＧ Ｗｅｉｃｈａｏ， ｅｔ ａｌ． Ｐａｒａｍｅｔｅｒ⁃ｉｎｄｕｃｅｄ ｓｔｏｃｈａｓｔｉｃ ｒｅｓｏｎａｎｃｅ ｉｎ ｍｕｌｔｉ⁃ｆｒｅｑｕｅｎｃｙ ｗｅａｋ ｓｉｇｎａｌ ｄｅｔｅｃｔｉｏｎ ｗｉｔｈ α ｓｔａｂｌｅ ｎｏｉｓｅ［Ｊ］． Ａｃｔａ Ｐｈｙｓｉｃａ Ｓｉｎｉｃａ，２０１３， ６２（２１）３４⁃４３． ［２０］ ＹＡＮＧ Ｘ Ｎ， ＹＡＮＧ Ｙ Ｆ， ＳＣＩＥＮＣＥ Ｓ Ｏ．Ｖｉｂｒａｔｉｏｎａｌ ｒｅｓｏｎａｎｃｅ ｉｎ ａｎ ａｓｙｍｍｅｔｒｉｃ ｂｉｓｔａｂｌｅ ｓｙｓｔｅｍ ｗｉｔｈ ｔｉｍｅ⁃ｄｅｌａｙ ｆｅｅｄｂａｃｋ ［ Ｊ ］．ＡｃｔａＰｈｙｓｉｃａＳｉｎｉｃａ， ２０１５， ６４ （ ７ ） ７０５０７⁃０７０５０７． 􀤖􀤖􀤖􀤖􀤖􀤖􀤖􀤖􀤖􀤖􀤖􀤖􀤖􀤖􀤖􀤖􀤖􀤖􀤖􀤖􀤖􀤖􀤖􀤖􀤖􀤖􀤖􀤖􀤖􀤖􀤖􀤖􀤖􀤖􀤖􀤖􀤖􀤖􀤖􀤖􀤖􀤖􀤖􀤖􀤖􀤖􀤖􀤖􀤖􀤖 （上接第 ５１ 页） ［ ３］ 徐明毅， 张勇传． 精细辛几何算法的误差估计［Ｊ］． 数学 物理学报， ２００６， ２６（２）３１４⁃３２０． ＸＵ Ｍｉｎｇｙｉ， ＺＨＡＮＧ Ｙｏｎｇｃｈｕａｎ．Ａｃｃｕｒａｃｙ ｅｓｔｉｍａｔｉｏｎ ｏｆ ｐｒｅｃｉｓｅ ｓｙｍｐｌｅｃｔｉｃ ｉｎｔｅｇｒａｔｉｏｎ ｍｅｔｈｏｄ［Ｊ］． Ａｃｔａ Ｍａｔｈｅｍａｔｉｃａ Ｓｃｉｅｎｔｉａ， ２００６， ２６（２）３１４⁃３２０． ［ ４］ ＭＯＬＥＲ Ｃ， ＶＡＮ Ｌ Ｃ． Ｎｉｎｅｔｅｅｎ ｄｕｂｉｏｕｓ ｗａｙｓ ｔｏ ｃｏｍｐｕｔｅ ｔｈｅ ｅｘｐｏｎｅｎｔｉａｌ ｏｆ ａ ｍａｔｒｉｘ， ｔｗｅｎｔｙ⁃ｆｉｖｅ ｙｅａｒｓ ｌａｔｅｒ［Ｊ］． ＳＩＡＭ Ｒｅｖｉｅｗ， ２００３， ４５（１） ３⁃４９． ［ ５］ 谭述君， 吴志刚， 钟万勰． 矩阵指数精细积分方法中参数 的自适应选择［Ｊ］． 力学学报， ２００９， ４１（６） ９６１⁃９６６． ＴＡＮ Ｓｈｕｊｕｎ， ＷＵ Ｚｈｉｇａｎｇ， ＺＨＯＮＧ Ｗａｎｘｉｅ．Ａｄａｐｔｉｖｅ ｓｅｌｅｃｔｉｏｎ ｏｆ ｐａｒａｍｅｔｅｒｓ ｆｏｒ ｐｒｅｃｉｓｅ ｃｏｍｐｕｔａｔｉｏｎ ｏｆ ｍａｔｒｉｘ ｅｘｐｏｎｅｎｔｉａｌ ｂａｓｅｄ ｏｎ Ｐａｄ ａｐｐｒｏｘｉｍａｔｉｏｎ ［ Ｊ ］．Ｃｈｉｎｅｓｅ Ｊｏｕｒｎａｌ ｏｆ Ｔｈｅｏｒｅｔｉｃａｌ ａｎｄ Ａｐｐｌｉｅｄ Ｍｅｃｈａｎｉｃｓ， ２００９， ４１（６） ９６１⁃９６６． ［ ６］ 向宇， 黄玉盈， 黄健强． 一种新型齐次扩容精细积分法 ［Ｊ］． 华中科技大学学报（自然科学版）， ２００２， ３０（１１） ７４⁃７６． ＸＩＡＮＧ Ｙｕ， ＨＵＡＮＧ Ｙｕｙｉｎｇ， ＨＵＡＮＧ Ｊｉａｎｑｉａｎｇ． Ａ ｍｅｔｈｏｄ ｏｆ ｈｏｍｏｇｅｎｉｚａｔｉｏｎ ｏｆ ｈｉｇｈ ｐｒｅｃｉｓｉｏｎ ｄｉｒｅｃｔ ｉｎｔｅｇｒａｔｉｏｎ ［ Ｊ］． Ｊｏｕｒｎａｌ ｏｆ Ｈｕａｚｈｏｎｇ Ｕｎｉｖｅｒｓｉｔｙ ｏｆ Ｓｃｉｅｎｃｅ ａｎｄ Ｔｅｃｈｎｏｌｏｇｙ （Ｎａｔｕｒａｌ Ｓｃｉｅｎｃｅ Ｅｄｉｔｉｏｎ）， ２００２， ３０（１１） ７４⁃７６． ［ ７］ 凌明祥， 韩宇航， 朱长春， 等． 非线性动力方程的三次样 条⁃增维精细算法［Ｊ］． 计算力学学报， ２０１４（６）７２９⁃７３４． ＬＩＮＧ Ｍｉｎｇｘｉａｎｇ， ＨＡＮ Ｙｕｈａｎｇ， ＺＨＵ Ｃｈａｎｇｃｈｕｎ， ｅｔ ａｌ． Ｉｎｃｒｅｍｅｎｔ⁃ｄｉｍｅｎｓｉｏｎａｌ ｐｒｅｃｉｓｅ ｉｎｔｅｇｒａｔｉｏｎ ｍｅｔｈｏｄ ｂａｓｅｄ ｏｎ ｃｕｂｉｃ ｓｐｌｉｎｅ ｉｎｔｅｒｐｏｌａｔｉｏｎ ｆｏｒ ｎｏｎｌｉｎｅａｒ ｄｙｎａｍｉｃ ｅｑｕａｔｉｏｎ［Ｊ］． Ｃｈｉｎｅｓｅ Ｊｏｕｒｎａｌ ｏｆ Ｃｏｍｐｕｔａｔｉｏｎａｌ Ｍｅｃｈａｎｉｃｓ， ２０１４ （６） ７２９⁃７３４． ［ ８］ 梅树立， 张森文， 徐加初． 结构非线性动力学方程的自适 应精细积分算法［Ｊ］． 华南理工大学学报（自然科学版）， ２００３， ３１（增刊 １）１３３⁃１３５． ＭＥＩ Ｓｈｕｌｉ， ＺＨＡＮＧ Ｓｅｎｗｅｎ， ＸＵ Ｊｉａｃｈｕ． Ａｄａｐｔｉｖｅ ｐｒｅｃｉｓｅ ｉｎｔｅｇｒａｔｉｏｎ ａｌｇｏｒｉｔｈｍ ｆｏｒ ｎｏｎｌｉｎｅａｒ ｄｙｎａｍｉｃａｌ ｅｑｕａｔｉｏｎ ［ Ｊ］． Ｊｏｕｒｎａｌ ｏｆ Ｓｏｕｔｈ Ｃｈｉｎａ Ｕｎｉｖｅｒｓｉｔｙ ｏｆ Ｔｅｃｈｎｏｌｏｇｙ （ Ｎａｔｕｒａｌ Ｓｃｉｅｎｃｅ Ｅｄｉｔｉｏｎ）， ２００３， ３１（Ｓｕｐ１）１３３⁃１３５． ［ ９］ ＭＡＴＨＥＷＳ Ｊ Ｈ， ＦＩＮＫ Ｋ Ｄ．Ｎｕｍｅｒｉｃａｌ ｍｅｔｈｏｄｓ ｕｓｉｎｇ ＭＡＴＬＡＢ［Ｍ］． Ｕｐｐｅｒ Ｓａｄｄｌｅ Ｒｉｖｅｒ， ＮＪ Ｐｅａｒｓｏｎ Ｐｒｅｎｔｉｃｅ Ｈａｌｌ， ２００４． ［１０］ ＬＩＮ Ｊ， ＳＨＥＮ Ｗ， ＷＩＬＬＩＡＭＳ Ｆ Ｗ． Ａｃｃｕｒａｔｅ ｈｉｇｈ⁃ｓｐｅｅｄ ｃｏｍｐｕｔａｔｉｏｎ ｏｆ ｎｏｎ⁃ｓｔａｔｉｏｎａｒｙ ｒａｎｄｏｍ ｓｔｒｕｃｔｕｒａｌ ｒｅｓｐｏｎｓｅ［Ｊ］． Ｅｎｇｉｎｅｅｒｉｎｇ Ｓｔｒｕｃｔｕｒｅｓ， １９９７， １９（７） ５８６⁃５９３． ０７振 动 与 冲 击 ２０２０ 年第 ３９ 卷