随机和区间非齐次线性哈密顿系统的比较研究及其应用 1)

doi:10.6052/0459-1879-19-348

随机和区间非齐次线性哈密顿系统的比较研究及其应用 ¹⁾

邱志平^,²⁾, 姜南

北京航空航天大学航空科学与工程学院, 北京 100191

COMPARATIVE STUDY OF STOCHASTIC AND INTERVAL NON-HOMOGENEOUS LINEAR HAMILTONIAN SYSTEMS AND THEIR APPLICATIONS ¹⁾

Qiu Zhiping^,²⁾, Jiang Nan

School of Aeronautic Science and Engineering, Beihang University, Beijing 100191, China

通讯作者: 2) 邱志平, 教授, 主要研究方向: 计算固体力学、结构动力学、哈密顿系统辛算法等. E-mail:zpqiu@buaa.edu.cn

收稿日期: 2019-12-9 接受日期: 2020-01-17 网络出版日期: 2020-01-18

基金资助:

1) 国家自然科学基金.  11772026
国防基础科研计划.  JCKY2016204B101
国防基础科研计划.  JCKY2018601B001

Received: 2019-12-9 Accepted: 2020-01-17 Online: 2020-01-18

作者简介 About authors

摘要

随着计算机技术的飞速发展,更高效、更稳定和长时间模拟能力更强的数值算法需求迫切.哈密顿系统辛算法与传统算法相比在稳定性和长期模拟方面具有显著优越性.但动力系统中不可避免地存在大量不同程度的不确定性,动力学分析中需要考虑这些不确定性的影响以确保合理有效性. 然而,目前考虑参数不确定性的哈密顿系统响应分析的研究基础还比较薄弱. 为此,本文考虑随机和区间参数不确定性,对两种不确定性非齐次线性哈密顿系统分析计算结果进行了比较研究,从而突破了传统哈密顿系统的局限性, 并应用于结构动力响应评估中. 首先,针对确定性非齐次线性哈密顿系统, 提出了考虑确定性扰动的参数摄动法;在此基础上, 分别提出了随机、区间非齐次线性哈密顿系统的参数摄动法,得到了它们响应界限的数学表达; 随后,用数学理论推导得到了区间响应范围包含随机响应范围的相容性结论; 最后,两个数值算例在较小时间步长下验证了所提方法在结构动力响应中的可行性和有效性,体现了随机、区间哈密顿系统响应结果之间的包络关系,并在较大时间步长下与传统方法相比较凸显了哈密顿系统辛算法的数值计算优势、与蒙特卡洛模拟方法相比较验证了所提方法的精度.

关键词： 哈密顿系统 ; 不确定性 ; 参数摄动法 ; 动力响应 ; 比较研究

Abstract

With the rapid development of computer technology, there is an urgent need for more efficient and more stable numerical algorithms with more powerful long-term simulation capabilities. Compared with the traditional algorithms, the symplectic algorithms of Hamiltonian systems have significant advantages in stability and long-term simulation. However, a variety of different degrees of uncertainties exist inevitably in the dynamic system, and the impacts of these uncertainties need to be considered in the dynamic analysis to ensure the rationality and effectiveness. Nevertheless, there has been very little research considering parameter uncertainties on the dynamic response analysis of Hamiltonian systems. For this reason, two kinds of uncertain non-homogeneous linear Hamiltonian systems are studied and compared in this paper breaking through the limitations of traditional Hamiltonian systems, where stochastic and interval parameter uncertainties are taken into account, and applied to the evaluation of structural dynamic response. Firstly, for the deterministic non-homogeneous linear Hamiltonian systems, a parameter perturbation method considering deterministic perturbations is proposed. On this basis, the parameter perturbation methods of stochastic and interval non-homogeneous linear Hamiltonian systems are proposed respectively, and the mathematical expressions of the bounds of their response are obtained. Then, the compatibility conclusion that the region of the dynamic response obtained by the interval method contains that obtained by the stochastic one is derived theoretically. Finally, two numerical examples verify the feasibility and effectiveness of the proposed method in structural dynamic response in a smaller time step, and reflect the envelope relationship between the numerical results of the response of stochastic and interval Hamiltonian systems. Also, in a larger time step, the numerical advantages of the symplectic algorithms of Hamiltonian systems are highlighted compared to the traditional algorithms and the accuracy of the proposed method is verified by comparison with the Monte Carlo simulation method.

Keywords： Hamiltonian system ; uncertainty ; parameter perturbation method ; dynamic response ; comparative study

PDF (1852KB) 元数据多维度评价相关文章导出 EndNote| Ris| Bibtex 收藏本文

本文引用格式

邱志平, 姜南. 随机和区间非齐次线性哈密顿系统的比较研究及其应用 ¹⁾. 力学学报[J], 2020, 52(1): 60-72 DOI:10.6052/0459-1879-19-348

Qiu Zhiping, Jiang Nan. COMPARATIVE STUDY OF STOCHASTIC AND INTERVAL NON-HOMOGENEOUS LINEAR HAMILTONIAN SYSTEMS AND THEIR APPLICATIONS ¹⁾. Chinese Journal of Theoretical and Applied Mechanics[J], 2020, 52(1): 60-72 DOI:10.6052/0459-1879-19-348

引言

随着计算机技术的飞速发展,研究人员越来越追求更高效、更稳定和长时间模拟能力更强的数值算法^[1].同一力学问题有牛顿力学、拉格朗日力学、哈密顿力学体系三种表示形式,其中一切可忽略耗散的物理过程都可以表示为某种哈密顿形式, 但由于求解途径不同,产生的计算结果可能是不等效的^[2].传统算法除少数例外, 都不是辛算法,不可避免地带有人为耗散性等歪曲体系特征的缺陷, 用于短期模拟尚可,用于长期跟踪则会导致结果严重失真^[3-4].而哈密顿系统辛算法却有保持体系结构的优点,在结构对称性和守恒性方面优于传统算法,特别在稳定性和长期跟踪能力上具有独特优越性^[2-5].

哈密顿系统辛算法的思想始于20世纪80年代Ruth^[6]和冯康^[7]的工作.1983年,Ruth^[6]构造了可分哈密顿系统的前三阶辛差分格式.1984年, 冯康^[7]首次系统地提出了哈密顿系统的辛算法,开创了哈密顿力学数值计算的新领域. 1988年,Sanz-Serna^[8]、Lasagni^[9]和Suris^[10]分别从不同角度给出了辛Runge-Kutta方法的判定条件.1993年,Sun^[11]对一般的哈密顿系统构造了辛分块Runge-Kutta方法.2000年前后,Bridges^[12]、Reich^[13]和Marsden等^[14]构建了多辛算法的理论框架.之后,变分积分子、离散梯度法、投影算法、分裂组合算法、对称算法等^[15]一系列方法不断提出并发展.近年来,无网格格式^[16]、连续级Runge-Kutta方法^[17]等辛算法也相继被提出.深入的理论分析和大量的数值实验令人信服地表明,辛算法在数值计算中具有显著优越性.

然而,动力系统中不可避免地存在大量的、不同程度的不确定性^[18].例如, 工程结构中的典型不确定性有: 建模过程中的简化操作导致所建模型存在误差,制造环境、材料多相特征等因素使弹性模量、泊松比等材料参数具有分散性,制造及安装误差使结构几何尺寸具有不确定性,测量条件、外部环境等因素使外载荷具有不确定性等^[19-20].这些不确定性将引起动力响应变化,响应不确定量有时甚至能达到参数不确定量的数倍^[21].此外, 不确定性的存在还会使系统的性质发生改变,保守系统若存在不确定性则不再保守. 根据产生机理不同,不确定性可分为随机不确定性和认知不确定性^[22-23]:随机不确定性是由自然变异和随机性而导致的不确定性, 以随机模型定量化;认知不确定性是指受知识水平和社会环境等因素制约而产生的认知上的不确定性,通常以区间模型定量化. 因此, 需要在哈密顿系统中考虑随机和区间不确定性的影响,以确保动力学分析计算的合理有效性.

1984年,Zambrini首先基于变分原理提出了Nelson随机力学在哈密顿力学下的运动方程,逐步建立了哈密顿力学下的随机力学体系. 2002年,Milstein等^[24]对一般的随机哈密顿系统构造了几类辛Runge-Kutta方法,开创了随机哈密顿系统辛算法这一全新的研究领域. 2007年,王丽瑾^[25]提出了随机变分积分子理论,使利用随机生成函数构造随机辛算法成为可能. 2009年,Bou-Rabee和Owhadi^[26]对随机哈密顿系统构造了随机变分积分子.近几年,丁效华课题组^[27-28]也对随机哈密顿系统进行了研究,并构造了随机辛Runge-Kutta、辛分块Runge-Kutta方法等. 此外,朱位秋^[29-30]研究了多自由度非线性随机动力学系统,在国际上首次提出了随机激励的耗散的哈密顿系统理论,构建了非线性随机动力学与控制的哈密顿理论框架.上述随机哈密顿系统相关研究主要聚焦随机白噪声激励,但没有考虑系统本身的参数随机性,并且考虑区间不确定性的哈密顿系统也未见有人研究.

哈密顿系统辛算法的工程结构应用方面,罗恩等^[31]建立了多自由度系统弹性动力学的相空间非传统哈密顿变分原理,提出了称之为辛时间子域法的辛算法, 精度和效率都具有明显优越性.邢誉峰课题组^[32-33]针对结构动力学方程构造了多种辛差分格式,得到了令人满意的结果.高强等^[34]用辛算法求解拉压刚度不同桁架的非线性动力问题,Li等^[35]构造了一类动力学初值问题的辛算法并应用于简谐振子和简支梁,Yang等^[36]用辛算法对超长细杆进行动力学数值模拟.也有学者对多辛算法的工程结构应用开展了相关研究,如梁^[37]、板^[38]、杆^[39]的动力响应等.然而, 用哈密顿系统求解含不确定性的结构动力响应还少有人研究.

本文针对随机和区间不确定性,对含参数不确定性的非齐次线性哈密顿系统的动力响应分析进行首次尝试,提出两种不确定性哈密顿系统的参数摄动法,突破传统哈密顿系统只适用于保守系统的局限性,并开展两种哈密顿系统不确定性响应的相容性研究,以期为结构动力响应评估提供更加有效稳定的数值计算方法.

1 确定性非齐次线性哈密顿系统

考虑哈密顿正则方程

(1) $ \begin{eqnarray} \label{eq1} \dfrac{\mbox{d}{z}}{\mbox{d}t}={J}^{-1}H_{z} \end{eqnarray} $

其中

(2) $ \begin{eqnarray} \label{eq2} && {z}=\left( {z_1 ,z_2 ,\cdots ,z_n ,z_{n+1} ,\cdots ,z_{2n} }\right)^{\rm T} \\ &&H_{z} =\left( {\dfrac{\partial H}{\partial z_1 },\dfrac{\partial H}{\partial z_2},\cdots ,\dfrac{\partial H}{\partial z_{2n} }} \right)^{\rm T},\quad {J}=\left( {{\begin{array}{c@{\ \ \ }c} {{\bf 0}} & {{I}_n } \\ {-{I}_n } & {{\bf 0}} \\ \end{array} }} \right) \end{eqnarray} $

其中, ${I}_n $是$n$阶单位矩阵; ${J}$是标准单位辛矩阵, 满足${J}^{-1}={J}^{\rm T}=-{J}$; $H$称为该系统的哈密顿函数.

哈密顿系统(1)称为线性的, 如果哈密顿函数$H$是${z}$的二次型

(3) $ \begin{eqnarray} \label{eq3} H\left( {z} \right)=\dfrac{1}{2}{z}^{\rm T}{Cz} \end{eqnarray} $

其中${C}$为对称矩阵, 即${C}^{\rm T}={C}$, 于是正则方程(1)能够表示为

(4) $ \begin{eqnarray} \label{eq4} \dfrac{\mbox{d}{z}}{\mbox{d}t}={Bz} \end{eqnarray} $

其中${B}={J}^{-1}{C}$是无穷小辛阵, 即满足${JB}+{B}^{\rm T}{J}={{\bf 0}}$.

非齐次线性哈密顿系统是在线性哈密顿系统基础上增加了非齐次项

(5) $ \begin{eqnarray} \label{eq5} \dfrac{\mbox{d}{z}}{\mbox{d}t}={Bz}+{F}(t) \end{eqnarray} $

其中${F}(t)$是与时间$t$有关的向量.

此时, 哈密顿函数${H}'$可以表示为

(6) $ \begin{eqnarray} \label{eq6} {H}'\left( {z} \right)=\dfrac{1}{2}{z}^{\rm T}{Cz}+{z}^{\rm T}{JF} \end{eqnarray} $

从而, 非齐次线性哈密顿系统(5)可以表示为哈密顿正则方程(1)的形式.

对于哈密顿系统, 其相空间配备着一个标准辛结构, 也就是一个闭的微分2形式

(7) $ \begin{eqnarray} \label{eq7} \omega =\sum \limits_{i=1}^n {\mbox{d}z_i \wedge \mbox{d}z_{n+i} } \end{eqnarray} $

哈密顿系统相流保持相空间的辛结构不变, 即

(8) $ \begin{eqnarray} \label{eq8} \dfrac{\mbox{d}\omega }{\mbox{d}t}=0 \end{eqnarray} $

2 含扰动非齐次线性哈密顿系统的参数摄动法

假设非齐次线性哈密顿系统(5)中的矩阵${B}$和向量${F}(t)$中元素与系统参数${b}=\left( {b_1 ,b_2 ,\cdots ,b_m }\right)^{\rm T}$有关, 此时, 哈密顿系统(5)可以写为

(9) $ \begin{eqnarray} \label{eq9} \dfrac{\mbox{d}{z}}{\mbox{d}t}={B}({b}){z}+{F}({b}) \end{eqnarray} $

当系统参数${b}=\left( {b_1 ,b_2 ,\cdots ,b_m } \right)^{\rm T}$存在扰动时, 矩阵${B}({b})$和向量${F}({b})$由标称值变化到扰动系统, 分析扰动对哈密顿系统响应${z}$的影响.

根据摄动理论, 引入小参数$\varepsilon $、参数${b}$、矩阵${B}$、向量${F}$和哈密顿系统的解${z}$可以分别写为摄动级数展开形式

(10) $ \begin{eqnarray} \label{eq10} \left.\begin{array}{l} {b}={b}_0 +\varepsilon {b}_1 +\varepsilon ^2{b}_2 +\cdots +\varepsilon ^p{b}_p +\cdots \\ {B}={B}_0 +\varepsilon {B}_1 +\varepsilon ^2{B}_2 +\cdots +\varepsilon ^p{B}_p +\cdots \\ {F}={F}_0 +\varepsilon {F}_1 +\varepsilon ^2{F}_2 +\cdots +\varepsilon ^p{F}_p +\cdots \\ {z}={z}_0 +\varepsilon {z}_1 +\varepsilon ^2{z}_2 +\cdots +\varepsilon ^p{z}_p +\cdots \\ \end{array}\right\} \end{eqnarray} $

其中, $\varepsilon $是小于单位1的小量, ${b}_0 $, ${B}_0 $, ${F}_0 $和${z}_0 $分别是${b}$, ${B}$, ${F}$和${z}$的标称值, ${b}_i $, ${B}_i $, ${F}_i $和${z}_i $ $(i=1,2,\cdots ,p)$分别是其第$i$阶摄动量. 含$\varepsilon$的项表示其与标称值相比是一个很小的量, 在这种情况下, 哈密顿系统的解${z}$只有小变化. 通过选择使$\left| \varepsilon \right|$足够小的参数可以使级数收敛.

将式(10)代入式(9), 得

(11) $ \begin{eqnarray} \label{eq11} &&{\mbox{d}\left( {{z}_0 +\varepsilon {z}_1 +\varepsilon ^2{z}_2 +\cdots +\varepsilon ^p{z}_p +\cdots } \right)} \Big/ {\mbox{d}t} =\\&&\qquad\left( {{B}_0 +\varepsilon {B}_1 +\varepsilon ^2{B}_2 +\cdots +\varepsilon ^p{B}_p +\cdots } \right) \cdot \\&&\qquad\left( {{z}_0 +\varepsilon {z}_1 +\varepsilon ^2{z}_2 +\cdots +\varepsilon ^p{z}_p +\cdots } \right) +\\&&\qquad\left( {{F}_0 +\varepsilon {F}_1 +\varepsilon ^2{F}_2 +\cdots +\varepsilon ^p{F}_p +\cdots } \right) \end{eqnarray} $

将式(11)展开, 比较$\varepsilon $的同次幂系数, 可得

(12) $ \begin{eqnarray} \label{eq12} \left.\begin{array}{l} \varepsilon ^0:\dfrac{\mbox{d}{z}_0 }{\mbox{d}t}={B}_0 {z}_0 +{F}_0 \\ \varepsilon ^1:\dfrac{\mbox{d}{z}_1 }{\mbox{d}t}={B}_0 {z}_1 +{B}_1 {z}_0 +{F}_1 \\ \varepsilon ^2:\dfrac{\mbox{d}{z}_2 }{\mbox{d}t}={B}_0 {z}_2 +{B}_1 {z}_1 +{B}_2 {z}_0 +{F}_2 \\ \cdots \cdots \\ \varepsilon ^n:\dfrac{\mbox{d}{z}_p }{\mbox{d}t}={B}_0 {z}_p +{B}_1 {z}_{p-1} +\cdots +{B}_{p-1} {z}_1 +\\ \qquad {B}_p {z}_0 +{F}_p \\ \end{array}\right\} \end{eqnarray} $

运用辛算法求解式(12)中第1式可以得到${z}$的标称部分, 即${z}_0 $. 本文中采用的辛算法均为欧拉中点格式. 对式(12)中其他式子的${B}_i $和${F}_i$ $(i=1,2,\cdots ,p)$, 可以通过Taylor展开获得, 即将${B}$和${F}$在${b}={b}_0 $处分别进行Taylor展开, 得到

(13a)$ \begin{eqnarray} && {B}\left( {b} \right)={B}\left( {{b}_0 +\varepsilon {b}_1 +\varepsilon ^2{b}_2 +\cdots +\varepsilon ^p{b}_p +\cdots } \right) =\\&&\qquad {B}_0 +\varepsilon \sum \limits_{j=1}^m {\dfrac{\partial {B}_0 }{\partial b_j }} b_{1j} +\varepsilon ^2\sum \limits_{j=1}^m {\dfrac{\partial {B}_0 }{\partial b_j }} b_{2j} +\cdots +\\&&\qquad\varepsilon ^p\sum \limits_{j=1}^m {\dfrac{\partial {B}_0 }{\partial b_j }} b_{pj} +\cdots \end{eqnarray} $

(13b)$ \begin{eqnarray} && {F}\left( {b} \right)={F}\left( {{b}_0 +\varepsilon {b}_1 +\varepsilon ^2{b}_2 +\cdots +\varepsilon ^p{b}_p +\cdots } \right)= \\&&\qquad {F}_0 +\varepsilon \sum \limits_{j=1}^m {\dfrac{\partial {F}_0 }{\partial b_j }} b_{1j} +\varepsilon ^2\sum \limits_{j=1}^m {\dfrac{\partial {F}_0 }{\partial b_j }} b_{2j} +\cdots +\\&&\qquad\varepsilon ^p\sum \limits_{j=1}^m {\dfrac{\partial {F}_0 }{\partial b_j }} b_{pj} +\cdots \end{eqnarray} $

其中$b_{ij}$ $(i=1,2,\cdots ,p;j=1,2,\cdots ,m)$是${b}_i $的分量.

由式(10)和式(13), 可知

(14) $ \begin{eqnarray} \label{eq14} {B}_i =\sum \limits_{j=1}^m {\dfrac{\partial {B}_0 }{\partial b_j }} b_{ij} ,\quad {F}_i =\sum \limits_{j=1}^m {\dfrac{\partial {F}_0 }{\partial b_j }} b_{ij} \end{eqnarray} $

将式(12)中第1式求得的${z}_0 $和式(14)求得的${B}_1 $, ${F}_1 $代入式(12)中第2式并利用辛算法求解, 可以求得${z}$的第1阶摄动量${z}_1 $; 进而可以通过依次求解式(12)中其他式子得到${z}_i$ $(i=2,3,\cdots ,p)$的值. 从而, ${z}$就可以按式(10)求得. 在每一步都采用辛算法求解保证计算结果能够保持体系结构特征, 避免传统算法带有的人为耗散性等歪曲体系特征的缺陷. 在实际计算中, 为了方便求解, 常常展开到第1阶摄动.

当参数${b}$存在的扰动为随机或区间不确定性时, 上述参数摄动法可以推广至求解随机或区间非齐次线性哈密顿系统, 详细求解过程如下面第3、4节所述.

3 随机非齐次线性哈密顿系统的参数摄动法

当系统参数${b}=\left( {b_1 ,b_2 ,\cdots ,b_m } \right)^{\rm T}$是随机变量时, 矩阵${B}$、向量${F}$和哈密顿系统的解${z}$也是随机的, 它们可以分别看作围绕确定性部分(均值)有一个随机小扰动. 因此, 基于前述摄动理论, 同样引入小参数$\varepsilon $, 将${b}$, ${B}$, ${F}$和${z}$分别表示为

(15) $ \begin{eqnarray} \label{eq15} \left.\begin{array}{l} {b}={b}_d +\varepsilon {b}_s ,\quad {B}={B}_d +\varepsilon {B}_s \\ {F}={F}_d +\varepsilon {F}_s ,\quad {z}={z}_d +\varepsilon {z}_s \\ \end{array}\right\} \end{eqnarray} $

其中, ${b}_d $, ${B}_d $, ${F}_d $和${z}_d $分别是${b}$, ${B}$, ${F}$和${z}$的确定性部分; ${b}_s $, ${B}_s $, ${F}_s $和${z}_s $分别是其随机部分, 且它们的均值均为0.

将式(15)代入式(9), 得

(16) $ \begin{eqnarray} \label{eq16} \dfrac{\mbox{d}\left( {{z}_d +\varepsilon {z}_s } \right)}{\mbox{d}t}=\left( {{B}_d +\varepsilon {B}_s } \right)\left( {{z}_d +\varepsilon {z}_s } \right)+\left( {{F}_d +\varepsilon {F}_s } \right) \end{eqnarray} $

将式(16)展开, 忽略$O\left( {\varepsilon ^2} \right)$高阶项, 并比较$\varepsilon $的同次幂系数, 可得

(17) $ \begin{eqnarray} \label{eq17} \left.\begin{array}{l} \varepsilon ^0:\dfrac{\mbox{d}{z}_d }{\mbox{d}t}={B}_d {z}_d +{F}_d\\ \varepsilon ^1:\dfrac{\mbox{d}{z}_s }{\mbox{d}t}={B}_d {z}_s +{B}_s {z}_d +{F}_s \\ \end{array}\right\} \end{eqnarray} $

利用辛算法求解式(17)中第1式可以求得${z}$的确定性部分${z}_d $, 即为${z}$的均值. 但无法由式(17)直接确定${z}$的随机部分${z}_s $, 需要进行变换后加以求解.

对式(15)求取数学期望, 有

(18) $ \begin{eqnarray} \label{eq18} \left.\begin{array}{l} E\left[ {b} \right]=E\left[ {{b}_d } \right]+E\left[ {\varepsilon {b}_s } \right]={b}_d \\ E\left[ {B} \right]=E\left[ {{B}_d } \right]+E\left[ {\varepsilon {B}_s } \right]={B}_d \\ E\left[ {F} \right]=E\left[ {{F}_d } \right]+E\left[ {\varepsilon {F}_s } \right]={F}_d \\ E\left[ {z} \right]=E\left[ {{z}_d } \right]+E\left[ {\varepsilon {z}_s } \right]={z}_d \\ \end{array}\right\} \end{eqnarray} $

由于${z}_d $是一个确定性的向量, 所以${z}_d $与${z}$, ${z}_d $与${z}_s $均相互独立, 从而可得

(19) $ \begin{eqnarray} \label{eq19} \left.\begin{array}{l} E\left[ {{zz}_d ^{\rm T}} \right]=E\left[ {z} \right]{z}_d ^{\rm T}={z}_d {z}_d ^{\rm T} \\ E\left[ {{z}_d {z}^{\rm T}} \right]={z}_d E\left[ {z} \right]^{\rm T}={z}_d {z}_d ^{\rm T} \\ E\left[ {{z}_d {z}_s ^{\rm T}} \right]={z}_d E\left[ {{z}_s } \right]^{\rm T}=0 \\ E\left[ {{z}_s {z}_d ^{\rm T}} \right]=E\left[ {{z}_s } \right]{z}_d ^{\rm T}=0 \\ \end{array}\right\} \end{eqnarray} $

因此, ${z}$的协方差矩阵可以写为

(20) $ \begin{eqnarray} \label{eq20} &&\mbox{Cov}\left[ {{z},{z}^{\rm T}} \right]=E\left[ {\left( {{z}-E\left[ {z} \right]} \right)\left( {{z}-E\left[ {z} \right]} \right)^{\rm T}} \right] =\\&&\qquad E\left[ {{zz}^{\rm T}-{z}E\left[ {z} \right]^{\rm T}-E\left[ {z} \right]{z}^{\rm T}+E\left[ {z} \right]E\left[ {z} \right]^{\rm T}} \right] =\\&&\qquad E\left[ {{zz}^{\rm T}-{zz}_d ^{\rm T}-{z}_d {z}^{\rm T}+{z}_d {z}_d ^{\rm T}} \right] =\\&&\qquad E\left[ {{zz}^{\rm T}} \right]-E\left[ {{zz}_d ^{\rm T}} \right]-E\left[ {{z}_d {z}} \right]^{\rm T}+E\left[ {{z}_d {z}_d ^{\rm T}} \right] =\\&&\qquad E\left[ {{zz}^{\rm T}} \right]-{z}_d {z}_d ^{\rm T} \end{eqnarray} $

式(20)中的$E\left[ {{zz}^{\rm T}} \right]$可以写为

(21) $ \begin{eqnarray} \label{eq21} && E\left[ {{zz}^{\rm T}} \right]=E\left[ {\left( {{z}_d +\varepsilon {z}_s } \right)\left( {{z}_d +\varepsilon {z}_s } \right)^{\rm T}} \right] =\\&&\qquad E\left[ {{z}_d {z}_d ^{\rm T}+\varepsilon {z}_d {z}_s ^{\rm T}+\varepsilon {z}_s {z}_d ^{\rm T}+\varepsilon ^2{z}_s {z}_s ^{\rm T}} \right] =\\&&\qquad E\left[ {{z}_d {z}_d ^{\rm T}} \right]+\varepsilon E\left[ {{z}_d {z}_s ^{\rm T}} \right]+\varepsilon E\left[ {{z}_s {z}_d ^{\rm T}} \right]+\varepsilon ^2E\left[ {{z}_s {z}_s ^{\rm T}} \right] =\\&&\qquad {z}_d {z}_d ^{\rm T}+\varepsilon ^2E\left[ {{z}_s {z}_s ^{\rm T}} \right] \end{eqnarray} $

将式(21)代入式(20), 可得协方差矩阵为

(22) $ \begin{eqnarray} \label{eq22} \mbox{Cov}\left[ {{z},{z}^{\rm T}} \right]=\varepsilon ^2E\left[ {{z}_s {z}_s ^{\rm T}} \right] \end{eqnarray} $

该协方差矩阵的对角元素表示各点的方差, 其他非对角元素表示各点间的协方差. 因此, ${z}$的各分量$z_i$ $(i=1,2,\cdots ,2n)$的方差为

(23) $ \begin{eqnarray} \label{eq23} D\left[ {z_i } \right]=\varepsilon ^2E\left[ {\left( {z_{si} } \right)^2} \right] \end{eqnarray} $

同理, ${b}$的各分量$b_j (j=1,2,\cdots ,m)$的方差为

(24) $ \begin{eqnarray} \label{eq24} D\left[ {b_j } \right]=\varepsilon ^2E\left[ {\left( {b_{sj} } \right)^2} \right] \end{eqnarray} $

将${z}$在${z}={z}_d $处进行一阶Taylor展开得到

(25) $ \begin{eqnarray} \label{eq25} {z}\left( {b} \right)={z}\left( {{b}_d +\varepsilon {b}_s } \right)={z}_d +\varepsilon \sum \limits_{j=1}^m {\dfrac{\partial {z}_d }{\partial b_j }} b_{sj} \end{eqnarray} $

其中, $b_{sj}$ $(j=1,2,\cdots ,m)$是${b}_s $的分量.

从而, ${z}$的随机部分${z}_s $可以表示为

(26) $ \begin{eqnarray} \label{eq26} {z}_s =\sum \limits_{j=1}^m {\dfrac{\partial {z}_d }{\partial b_j }} b_{sj} \end{eqnarray} $

同理, ${B}_s $和${F}_s $可以表示为

(27) $ \begin{eqnarray} \label{eq27} {B}_s =\sum \limits_{j=1}^m {\dfrac{\partial {B}_d }{\partial b_j }} b_{sj} ,\quad {F}_s =\sum \limits_{j=1}^m {\dfrac{\partial {F}_d }{\partial b_j }} b_{sj} \end{eqnarray} $

将式(26) ~式(27)代入式(17)中第2式, 得

(28) $ \begin{eqnarray} \label{eq28} &&\dfrac{\mbox{d}}{\mbox{d}t}\left( {\sum \limits_{j=1}^m {\dfrac{\partial {z}_d }{\partial b_j }} b_{sj} } \right) ={B}_d \left( {\sum \limits_{j=1}^m {\dfrac{\partial {z}_d }{\partial b_j }} b_{sj} } \right)+\\&&\qquad\left( {\sum \limits_{j=1}^m {\dfrac{\partial {B}_d }{\partial b_j }} b_{sj} } \right){z}_d +\sum \limits_{j=1}^m {\dfrac{\partial {F}_d }{\partial b_j }} b_{sj} \end{eqnarray} $

运用辛算法求解式(28)可以得到$\left( {\sum \limits_{j=1}^m {\dfrac{\partial {z}_d }{\partial b_j }} b_{sj} } \right)$的值, 即为${z}_s $, 将其分量$z_{si} =\sum \limits_{j=1}^m {\dfrac{\partial z_{di} }{\partial b_j }} b_{sj} $代入式(23), 可以得到

(29) $ \begin{eqnarray} \label{eq29} &&D\left[ {z_i } \right]=\varepsilon ^2E\left[ {\left( {z_{si} } \right)^2} \right]=\varepsilon ^2E\left[ {\left( {\sum \limits_{j=1}^m {\dfrac{\partial z_{di} }{\partial b_j }} b_{sj} } \right)^2} \right] =\\&&\qquad \varepsilon ^2\sum \limits_{j=1}^m {\left( {\dfrac{\partial z_{di} }{\partial b_j }} \right)^2} E\left[ {\left( {b_{sj} } \right)^2} \right] +\\&&\qquad\varepsilon ^2\sum \limits_{k=1}^m {\sum \limits_{l=1}^m {\dfrac{\partial z_{di} }{\partial b_k }} \dfrac{\partial z_{di} }{\partial b_l }\mbox{Cov}\left( {b_k ,b_l } \right)} \end{eqnarray} $

其中$\mbox{Cov}\left( {b_k ,b_l } \right)$表示$b_k $和$b_l $的协方差.

若参数$b_j $是相互独立的, $\mbox{Cov}\left( {b_k ,b_l } \right)$为0, 则式(29)可以简化为

(30) $ \begin{eqnarray} \label{eq30} D\left[ {z_i } \right]=\varepsilon ^2\sum \limits_{j=1}^m {\left( {\dfrac{\partial z_{di} }{\partial b_j }} \right)^2} E\left[ {\left( {b_{sj} } \right)^2} \right] \end{eqnarray} $

将式(24)代入式(30), 就可以得到${z}$的各分量$z_i (i=1,2,\cdots ,2n)$的方差

(31) $ \begin{eqnarray} \label{eq31} D\left[ {z_i } \right]=\sum \limits_{j=1}^m {\left( {\dfrac{\partial z_{di} }{\partial b_j }} \right)^2} D\left[ {b_j } \right] \end{eqnarray} $

因此, $z_i$ $(i=1,2,\cdots ,2n)$的标准差为

(32) $ \begin{eqnarray} \label{eq32} && \sigma \left[ {z_i } \right]=\sqrt {D\left[ {z_i } \right]} =\sqrt {\sum \limits_{j=1}^m {\left( {\dfrac{\partial z_{di} }{\partial b_j }} \right)^2} D\left[ {b_j } \right]} =\\&&\qquad \sqrt {\sum \limits_{j=1}^m {\left( {\dfrac{\partial z_{di} }{\partial b_j }\sigma \left[ {b_j } \right]} \right)^2} } \end{eqnarray} $

其中$\sigma \left[ {b_j } \right]$是$b_j$ $(j=1,2,\cdots ,m)$的标准差.

在计算求解过程中, 求解式(17)中第1式求得${z}$的均值${z}_d $和求解式(28)求得${z}$的随机部分${z}_s $都采用了辛算法,确保计算结果能够保持体系结构特征.

令$k$为正整数, 则对$z_i$ $(i=1,2,\cdots ,2n)$而言, 距其均值的距离为$k$倍标准差的范围为

(33) $ \begin{eqnarray} \label{eq33} z_i^I =\left[ \underline{v}_i ,\bar{{v}}_i \right]=\left[ {z_{di} -k\sigma \left[ {z_i } \right],z_{di} +k\sigma \left[ {z_i } \right]} \right] \end{eqnarray} $

其中范围的上界和下界分别为

(34) $ \begin{eqnarray} \label{eq34} \left.\begin{array}{l} \bar{{v}}_i =z_{di} +k\sigma \left[ {z_i } \right]=z_{di} +k\sqrt {\sum \limits_{j=1}^m {\left( {\dfrac{\partial z_{di} }{\partial b_j }\sigma \left[ {b_j } \right]} \right)^2}} \\ \underline{v}_i =z_{di} -k\sigma \left[ {z_i } \right]=z_{di} -k\sqrt {\sum \limits_{j=1}^m {\left( {\dfrac{\partial z_{di} }{\partial b_j }\sigma \left[ {b_j } \right]} \right)^2} } \\ \end{array}\ \ \right\} \end{eqnarray} $

4 区间非齐次线性哈密顿系统的参数摄动法

当系统参数${b}=\left( {b_1 ,b_2 ,\cdots ,b_m } \right)^{\rm T}$是区间变量时, 即参数${b}$在一个区间向量内取值,矩阵${B}$、向量${F}$和哈密顿系统的解${z}$也分别在一个区间范围内取值

(35) $ \begin{eqnarray} \label{eq35} \left.\begin{array}{l} {b}\in {b}^I=\left[ {\underline{{b}},{\bar{{b}}}} \right],\quad {B}\in {B}^I=\left[ {\underline{{B}},{\bar{{B}}}} \right] \\ {F}\in {F}^I=\left[ {\underline{{F}},{\bar{{F}}}} \right],\quad {z}\in {z}^I=\left[ {\underline{{z}},{\bar{{z}}}} \right] \\ \end{array}\right\} \end{eqnarray} $

其中, ${b}^I$, ${F}^I$和${z}^I$是区间向量, ${B}^I$是区间矩阵, $\underline{{b}}$, $\underline{{B}}$, $\underline{{F}}$和$\underline{{z}}$分别是${b}$, ${B}$, ${F}$和${z}$的下界, ${\bar{{b}}}$, ${\bar{{B}}}$, ${\bar{{F}}}$和${\bar{{z}}}$分别是其上界.

此时, 哈密顿方程(9)可写为

(36) $ \begin{eqnarray} \label{eq36} \dfrac{\mbox{d}{z}^I\left( {{b}^I} \right)}{\mbox{d}t}={B}^I\left( {{b}^I} \right){z}^I\left( {{b}^I} \right)+{F}^I\left( {{b}^I} \right) \end{eqnarray} $

式(36)是区间非齐次线性哈密顿方程.

${b}^I$, ${B}^I$, ${F}^I$和${z}^I$的中值和半径分别为

(37) $ \begin{eqnarray} \label{eq37} \left.\begin{array}{l} {b}^c={\left( {{\bar{{b}}}+\underline{{b}}} \right)} / 2,\quad \Delta {b}={\left( {{\bar{{b}}}-\underline{{b}}} \right)} \Big/ 2 \\ {B}^c={\left( {{\bar{{B}}}+\underline{{B}}} \right)} / 2,\quad \Delta {B}={\left( {{\bar{{B}}}-\underline{{B}}} \right)} \Big/ 2 \\ {F}^c={\left( {{\bar{{F}}}+\underline{{F}}} \right)} / 2,\quad \Delta {F}={\left( {{\bar{{F}}}-\underline{{F}}} \right)} \Big/ 2\\ {z}^c={\left( {{\bar{{z}}}+\underline{{z}}} \right)} / 2,\quad \Delta {z}={\left( {{\bar{{z}}}-\underline{{z}}} \right)} \Big/ 2 \\ \end{array}\right\} \end{eqnarray} $

利用区间中心表示法, ${b}^I$, ${B}^I$, ${F}^I$和${z}^I$可以分别表示为

(38) $ \begin{eqnarray} \label{eq38} \left.\begin{array}{l} {b}^I={b}^c+\Delta {b}^I,\quad {B}^I={B}^c+\Delta {B}^I \\ {F}^I={F}^c+\Delta {F}^I,\quad {z}^I={z}^c+\Delta {z}^I \\ \end{array}\right\} \end{eqnarray} $

其中

(39) $ \begin{eqnarray} \label{eq39} \left.\begin{array}{l} \Delta {b}^I=\left[ {-\Delta {b},\Delta {b}} \right],\quad \Delta {B}^I=\left[ {-\Delta {B},\Delta {B}} \right] \\ \Delta {F}^I=\left[ {-\Delta {F},\Delta {F}} \right],\quad \Delta {z}^I=\left[ {-\Delta {z},\Delta {z}} \right] \\ \end{array}\right\} \end{eqnarray} $

把式(38)代入式(36), 得

(40) $ \begin{eqnarray} \label{eq40} \dfrac{\mbox{d}}{\mbox{d}t}\left( {{z}^c+\Delta {z}^I} \right)=\left( {{B}^c+\Delta {B}^I} \right)\left( {{z}^c+\Delta {z}^I} \right)+ {{F}^c+\Delta {F}^I} \end{eqnarray} $

如果将$\Delta {b}^I$, $\Delta {B}^I$, $\Delta {F}^I$和$\Delta {z}^I$分别看作围绕${b}^c$, ${B}^c$, ${F}^c$和${z}^c$的扰动, 则可以采用第2节所述参数摄动法求解区间哈密顿方程(40).

按照区间的含义, 引入小参数$\varepsilon $, 式(40)可以表示为: 由于参数${b}$存在小扰动$\delta {b}$, 导致${B}$, ${F}$和${z}$产生小扰动$\delta {B}$, $\delta {F}$和$\delta {z}$, 且满足

(41) $ \begin{eqnarray} \label{eq41} \left.\begin{array}{l} -\Delta {b}\leqslant \delta {b}\leqslant \Delta {b},\quad -\Delta {B}\leqslant \delta {B}\leqslant \Delta {B} \\ -\Delta {F}\leqslant \delta {F}\leqslant \Delta {F},\quad -\Delta {z}\leqslant \delta {z}\leqslant \Delta {z} \\ \end{array}\right\} \end{eqnarray} $

条件下的扰动方程的形式

(42) $ \begin{eqnarray} \label{eq42} \dfrac{\mbox{d}}{\mbox{d}t}\left( {{z}^c+\varepsilon \delta {z}} \right)=\left( {{B}^c+\varepsilon \delta {B}} \right)\left( {{z}^c+\varepsilon \delta {z}} \right)+ {{F}^c+\varepsilon \delta {F}} \end{eqnarray} $

式(41)和式(42)所表示的问题可以理解为: 在参数中值${b}^c$已知, 从而能够确定中值${B}^c$和${F}^c$, 而小扰动$\delta {b}$的具体取值未知但其取值范围(41)已知, 小扰动$\delta {B}$和$\delta {F}$的具体取值也未知但其取值范围(41)可以确定的情况下, 确定哈密顿系统的解${z}$的界限.

展开式(42)并比较$\varepsilon $的同次幂系数, 可得

(43) $ \begin{eqnarray} \label{eq43} \left.\begin{array}{l} \varepsilon ^0:\dfrac{\mbox{d}{z}^c}{\mbox{d}t}={B}^c{z}^c+{F}^c \\ \varepsilon ^1:\dfrac{\mbox{d}\delta {z}}{\mbox{d}t}={B}^c\delta {z}+\delta {Bz}^c+\delta {F} \\ \end{array}\right\} \end{eqnarray} $

运用辛算法求解式(43)中第1式可以求得${z}$的中值${z}^c$. 由区间扩张, 式(43)中第2式可写为

(44) $ \begin{eqnarray} \label{eq44} \dfrac{\mbox{d}\Delta {z}}{\mbox{d}t}={B}^c\Delta {z}+\Delta {Bz}^c+\Delta {F} \end{eqnarray} $

将${z}$在${z}={z}^c$处进行一阶Taylor展开得到

(45) $ \begin{eqnarray} \label{eq45} {z}\left( {b} \right)={z}\left( {{b}^c+\varepsilon \delta {b}} \right)={z}^c+\varepsilon \sum \limits_{j=1}^m {\dfrac{\partial {z}^c}{\partial b_j }} \delta b_j \end{eqnarray} $

其中, $\delta b_j$ $(j=1,2,\cdots ,m)$是$\delta {b}$的分量.

从而可得$\delta {z}$的表达式

(46) $ \begin{eqnarray} \label{eq46} \delta {z}=\sum \limits_{j=1}^m {\dfrac{\partial {z}^c}{\partial b_j }} \delta b_j \end{eqnarray} $

同理可得

(47) $ \begin{eqnarray} \label{eq47} \delta {B}=\sum \limits_{j=1}^m {\dfrac{\partial {B}^c}{\partial b_j }} \delta b_j ,\quad \delta {F}=\sum \limits_{j=1}^m {\dfrac{\partial {F}^c}{\partial b_j }\delta } b_j \end{eqnarray} $

式(46)和式(47)的区间扩张形式为

(48) $ \begin{eqnarray} \label{eq48} \left.\begin{array}{l} \Delta {z}=\sum \limits_{j=1}^m {\left| {\dfrac{\partial {z}^c}{\partial b_j }} \right|} \Delta b_j ,\quad \Delta {B}=\sum \limits_{j=1}^m {\left| {\dfrac{\partial {B}^c}{\partial b_j }} \right|} \Delta b_j \\ \Delta {F}=\sum \limits_{j=1}^m {\left| {\dfrac{\partial {F}^c}{\partial b_j }} \right|\Delta } b_j \\ \end{array}\right\} \end{eqnarray} $

其中, $\Delta b_j$ $(j=1,2,\cdots ,m)$是$\Delta {b}$的分量.

将式(48)代入式(44), 得

(49) $ \begin{eqnarray} \label{eq49} && \dfrac{\mbox{d}}{\mbox{d}t}\left( {\sum \limits_{j=1}^m {\left| {\dfrac{\partial {z}^c}{\partial b_j }} \right|} \Delta b_j } \right) ={B}^c\left( {\sum \limits_{j=1}^m {\left| {\dfrac{\partial {z}^c}{\partial b_j }} \right|} \Delta b_j } \right)+\\&&\qquad\left( {\sum \limits_{j=1}^m {\left| {\dfrac{\partial {B}^c}{\partial b_j }} \right|} \Delta b_j } \right){z}^c+ \sum \limits_{j=1}^m {\left| {\dfrac{\partial {F}^c}{\partial b_j }} \right|\Delta } b_j \end{eqnarray} $

利用辛算法求解式(49)可以得到$\sum \limits_{j=1}^m {\left| {\dfrac{\partial {z}^c}{\partial b_j }} \right|} \Delta b_j $的值, 即为$\Delta {z}$, 从而得到

(50) $ \begin{eqnarray} \label{eq50} {z}^I=\left[ {\underline{{z}},{\bar{{z}}}} \right]={z}^c+\Delta {z}^I={z}^c+\left[ {-\Delta {z},\Delta {z}} \right] \end{eqnarray} $

其中

(51) $ \begin{eqnarray} \label{eq51} \Delta {z}^I=\left[ {-\Delta {z},\Delta {z}} \right]=\left[ {-\sum \limits_{j=1}^m {\left| {\dfrac{\partial {z}^c}{\partial b_j }} \right|} \Delta b_j ,\sum \limits_{j=1}^m {\left| {\dfrac{\partial {z}^c}{\partial b_j }} \right|} \Delta b_j } \right] \end{eqnarray} $

在计算求解过程中, 求解式(43)中第1式求得${z}$的中值${z}^c$和求解式(49)求得${z}$的半径$\Delta {z}$都采用了辛算法,同样确保计算结果能够保持体系结构特征.

由区间相等的定义可得${z}$的上界和下界分别为

(52) $ \begin{eqnarray} \label{eq52} \left.\begin{array}{l} {\bar{{z}}}={z}^c+\Delta {z}={z}^c+\sum \limits_{j=1}^m {\left| {\dfrac{\partial {z}^c}{\partial b_j }} \right|} \Delta b_j\\ \underline{{z}}={z}^c-\Delta {z}={z}^c-\sum \limits_{j=1}^m {\left| {\dfrac{\partial {z}^c}{\partial b_j }} \right|} \Delta b_j \\ \end{array}\right\} \end{eqnarray} $

${z}$的分量形式$z_i$ $(i=1,2,\cdots ,2n)$的上界和下界分别为

(53) $ \begin{eqnarray} \label{eq53} \left.\begin{array}{l} \bar{{z}}_i =z_i^c +\Delta z_i =z_i^c +\sum \limits_{j=1}^m {\left| {\dfrac{\partial z_i^c }{\partial b_j }} \right|} \Delta b_j \\ \underline{z}_i =z_i^c -\Delta z_i =z_i^c -\sum \limits_{j=1}^m {\left| {\dfrac{\partial z_i^c }{\partial b_j }} \right|} \Delta b_j \\ \end{array}\right\} \end{eqnarray} $

5 随机和区间非齐次线性哈密顿系统比较

不确定性是客观存在的, 无论是随机方法还是区间方法, 只是描述不确定性的不同形式, 不能从本质上改变不确定性对响应的影响规律, 利用两种方法得到的不确定性响应理应具有相容性. 因此, 本节对随机和区间非齐次线性哈密顿系统的分析结果进行比较, 探究两者响应界限的包含关系.

假定参数${b}$的区间范围由概率统计信息获取, 即可以表示为距其均值距离为$k$倍标准差的形式

(54) $ \begin{eqnarray} \label{eq54} {b}^I =\left[ {\underline{{b}},{\bar{{b}}}} \right]=\left[ {{b}_d -k\sigma \left[ {b} \right],{b}_d +k\sigma \left[ {b} \right]} \right] \end{eqnarray} $

其中, $k$为正整数, $\sigma \left[ {b} \right]$是${b}$的标准差, ${b}$的上界和下界分别为

(55) $ \begin{eqnarray} \label{eq55} {\bar{{b}}}={b}_d +k\sigma \left[ {b} \right],\quad \underline{{b}}={b}_d -k\sigma \left[ {b} \right] \end{eqnarray} $

式(54)的分量形式$b_j^I$ $(j=1,2,\cdots ,m)$为

(56) $ \begin{eqnarray} \label{eq56} b_j^I =\left[ {\underline{b}_j ,\bar{{b}}_j } \right]=\left[ {b_{dj} -k\sigma \left[ {b_j } \right],b_{dj} +k\sigma \left[ {b_j } \right]} \right] \end{eqnarray} $

由式(54) ~式(55), 可以得到参数${b}$的区间中值和半径与随机确定性部分和标准差之间存在关系

(57) $ \begin{eqnarray} \label{eq57} {b}^c={b}_d ,\ b_j^c =b_{dj} ,\ \Delta {b}=k\sigma \left[ {b} \right],\ \Delta b_j =k\sigma \left[ {b_j } \right] \end{eqnarray} $

哈密顿系统的解${z}$的区间中值与随机确定性部分同样存在关系

(58) $ \begin{eqnarray} \label{eq58} \left.\begin{array}{l} {z}^c={z}\left( {{b}^c} \right)={z}\left( {{b}_d } \right)={z}_d \\ z_i^c =z_i \left( {b_j^c } \right)=z_i \left( {b_{dj} } \right)=z_{di}, \quad i=1,2,\cdots ,2n \\ \end{array}\right\} \end{eqnarray} $

将式(57)的分量形式代入式(53), 得到$z_i$ $(i=1,2,\cdots ,2n)$的上界和下界

(59) $ \bar{{z}}_i =z_{di}\!+\!\sum \limits_{j=1}^m {\left| {\dfrac{\partial z_{di} }{\partial b_j }} \right|} k\sigma \left[ {b_j } \right], \underline{z}_i =z_{di} \!-\!\sum \limits_{j=1}^m {\left| {\dfrac{\partial z_{di} }{\partial b_j }} \right|} k\sigma \left[ {b_j } \right] $

对于参数$b_j $, 切比雪夫不等式成立

(60) $ \begin{eqnarray} \label{eq60} \sum \limits_{j=1}^m {b_j } \geqslant \sqrt {\sum \limits_{j=1}^m {b_j^2 } } \end{eqnarray} $

基于不等式(60), 对于表达式$\left| {\dfrac{\partial z_{di}}{\partial b_j }} \right|k\sigma \left[ {b_j } \right]$, 有

(61) $ \begin{eqnarray} \label{eq61} &&\sum \limits_{j=1}^m {\left| {\dfrac{\partial z_{di} }{\partial b_j }} \right|} k\sigma \left[ {b_j } \right]\geqslant \sqrt {\sum \limits_{j=1}^m {\left( {\left| {\dfrac{\partial z_{di} }{\partial b_j }} \right|k\sigma \left[ {b_j } \right]} \right)^2} } =\\&&\qquad k\sqrt {\sum \limits_{j=1}^m {\left( {\dfrac{\partial z_{di} }{\partial b_j }\sigma \left[ {b_j } \right]} \right)^2} } \end{eqnarray} $

由不等式(61), 可以得到由随机方法确定的上下界(34)和由区间方法确定的上下界(59)存在关系

(62) $ \begin{eqnarray} \label{eq62} \underline{z}_i \leqslant \underline{v}_i \leqslant \bar{{v}}_i \leqslant \bar{{z}}_i, \quad i=1,2,\cdots ,2n \end{eqnarray} $

式(62)表示, 在由概率统计信息确定不确定性参数的区间范围的情况下, 对于不确定性非齐次线性哈密顿系统, 由区间方法获得的哈密顿系统的解的范围比由随机方法获得的范围大, 即区间方法得到的上界比随机方法得到的上界大, 而区间方法得到的下界比随机方法得到的下界小.

6 数值算例

为了验证所提方法在结构动力响应中的可行性和有效性, 本节给出两个数值算例, 包括悬臂梁和复合材料层合板, 并将本文所提随机、区间方法(分别简记为SHPM、IHPM)计算结果与传统随机、区间方法(分别简记为TSM、TIM)计算结果相比较.

6.1 悬臂梁

考虑如图1所示11节点、10单元悬臂梁在正弦激励作用下的动力响应. 梁长$L=1$ m, 横截面积$A=2$ cm$^2$, 横截面的惯性矩$I_z =2$ cm$^4$, 材料泊松比$\nu =0.3$. 正弦激励$P(t)=-p\sin (1600\pi t)$ N作用在节点3的竖直方向, 初始条件为${\dot{{x}}}(0)=0,{x}(0)=0$.

图 1

新窗口打开| 下载原图ZIP| 生成PPT

图 1 11节点、10单元悬臂梁

Fig. 1 A cantilever beam with 11 nodes, 10 elements

由于材料中不可避免的分散性及测量误差, 材料的弹性模量、密度和正弦激励幅值均具有不确定性, 假设它们所在的区间范围$I$如表1所示; 同时假设它们在区间范围内服从正态分布, 均值$\mu $和标准差$\sigma $也如表1所示.

表 1 材料弹性模量、密度和正弦激励幅值的不确定性

Table 1 Uncertainties of elastic modulus, density of the material and the amplitude of the harmonic excitation

新窗口打开| 下载CSV

悬臂梁整体运动微分方程为

(63) $ \begin{eqnarray} \label{eq63} {M\ddot{{x}}}+{Kx}={F}\left( t \right) \end{eqnarray}$

其中, ${M}$是悬臂梁整体质量矩阵, 与材料密度$\rho $有关; ${K}$是整体刚度矩阵, 与材料弹性模量$E$有关; ${F}(t)$是整体载荷向量, 与正弦激励$P(t)$有关.

令${y}={M\dot{{x}}}$, 则整体运动方程(63)可以表示为

(64) $ \begin{eqnarray} \label{eq64} \left( {{\begin{array}{*{20}c} {{\dot{{x}}}} \\ {{\dot{{y}}}} \\ \end{array} }} \right)=\left( {{\begin{array}{*{20}c} {{\bf 0}} & {{M}^{-1}} \\ {-{K}} & {{\bf 0}} \\ \end{array} }} \right)\left( {{\begin{array}{*{20}c} {x} \\ {y} \\ \end{array} }} \right)+\left( {{\begin{array}{*{20}c} {{\bf 0}} \\ {{F}\left( t \right)} \\ \end{array} }} \right) \end{eqnarray}$

由于${B}=\left(\begin{array}{cc} {\bf 0} & {M}^{-1} \\ -{K} & {\bf 0} \\ \end{array}\right)$满足${JB}+{B}^{\rm T}{J}={{\bf 0}}$, ${B}$即为无穷小辛阵, 式(64)所表示系统为非齐次线性哈密顿系统.

首先考察在时间步长$\Delta t=40$ $\mu$s下利用辛算法求解哈密顿系统(64)和直接求解方程(63)方法计算节点11的响应标称值.利用辛算法求解哈密顿系统(64)计算得到的$t=0\sim 0.1$ s内的响应标称值曲线如图2(a)所示, 整体呈周期性变化, 但直接求解方程(63)得到的响应标称值很短时间内即发散, 如图2(b)所示, 只有当时间步长足够小, 如$\Delta t=2$ $\mu$s时,直接求解方程(63)才能得到和利用辛算法求解哈密顿系统(64)相同的结果. 这一稳定性差异反映了哈密顿系统辛算法能够保持体系结构特征,体现出利用哈密顿系统辛算法求解微分方程的优越性.

图 2

新窗口打开| 下载原图ZIP| 生成PPT

图 2 时间步长$\Delta t=40$ $\mu$s下利用不同算法计算得到的节点11的响应标称值

Fig. 2 The nominal value of the response of node 11 obtained by different algorithms with time step $\Delta t=40$ $\mu$s

图 2

新窗口打开| 下载原图ZIP| 生成PPT

图 2 时间步长$\Delta t=40$ $\mu$s下利用不同算法计算得到的节点11的响应标称值(续)

Fig. 2 The nominal value of the response of node 11 obtained by different algorithms with time step $\Delta t=40$ $\mu$s (continued)

在同一时间步长$\Delta t=2$ $\mu$s下本文所提随机、区间非齐次线性哈密顿系统的参数摄动法和传统随机、区间方法计算得到的节点11在$t=0\sim 5.0$ ms内的响应曲线如图3所示, 其中图3(a) ~图3(b)分别为两种随机、区间方法得到的响应曲线, 图3(c)为本文所提随机、区间方法计算得到的响应曲线, 响应标称值也绘制于图3中.

图 3

新窗口打开| 下载原图ZIP| 生成PPT

图 3 本文所提方法和传统方法计算得到的节点11的响应曲线

Fig. 3 The response curve of node 11 obtained by SHPM, IHPM, TSM and TIM

本文所提随机、区间方法和传统随机、区间方法计算得到的节点11在$t=3.5$ ms时刻的位移上下界如表2所示.

表 2 本文所提随机、区间方法和传统随机、区间方法计算得到的节点11在$t=3.5$ ms时刻的位移上下界

Table 2 The upper and lower bounds on the displacement of node 11 at $t=3.5$ ms obtained by SHPM, IHPM TSM and TIM

新窗口打开| 下载CSV

由图3和表2可知, 本文所提随机、区间方法得到的响应上下界曲线分别与传统随机、区间方法所得响应上下界曲线均几乎完全重合,结果非常相近, 验证了所提随机、区间方法的准确性和有效性. 此外, 本文所提区间方法得到的响应区间范围包含本文所提随机方法得到的响应区间范围, 即区间方法所得响应上界大于随机方法所得响应上界,而区间方法下界小于随机方法下界, 这一现象与前述理论推导相符.

前述内容验证了在较大时间步长$\Delta t=40$ $\mu$s下哈密顿系统辛算法相较于直接求解运动微分方程方法计算响应标称值的有效性和优越性, 也验证了在较小时间步长$\Delta t=2$ $\mu$s下本文所提随机、区间方法所得响应区间范围的准确性. 为了进一步检验本文所提随机、区间方法在较大时间步长下所得响应区间范围也能保持较高精度, 下面将本文所提随机、区间方法在时间步长$\Delta t=40$ $\mu$s下得到的响应区间范围与蒙特卡洛模拟方法(简记为MCM)得到的响应区间范围进行比较. 其中, 蒙特卡洛模拟在较小时间步长$\Delta t=2$ $\mu$s下进行, 样本点数目设置为$1.0\times10^{5}$, 对于每一样本点都采用辛算法求解,从而可将蒙特卡洛模拟结果视为精确值. 本文所提随机、区间方法和蒙特卡洛模拟方法计算得到的节点11在$t=0\sim 5.0$ ms内的响应曲线如图4所示.

图 4

新窗口打开| 下载原图ZIP| 生成PPT

图 4 本文所提随机、区间方法和蒙特卡洛模拟方法计算得到的节点11的响应曲线

Fig. 4 The response curve of node 11 obtained by SHPM, IHPM and MCM

由图4所示, 本文所提随机方法与蒙特卡洛模拟方法得到的响应上下界曲线差别微小, 尽管由于较大时间步长原因导致精度有一定下降, 但仍得到了较为满意的结果. 这一现象体现出本文所提方法在较大时间步长下也能得到具有较高精度的结果, 再次验证了本文所提方法的优势. 而对于本文所提区间方法, 由于区间扩张效应, 得到的响应区间范围包含本文所提随机方法和蒙特卡洛模拟方法得到的响应区间范围,与前述理论推导相一致.

6.2 复合材料层合板

考虑如图5所示的边长为$L=1.0$的四边固支的复合材料层合板结构, 其由5层正交各向异性材料铺设而成, 铺设角度为 (0$^\circ$, 90$^\circ$, 0$^\circ$, 90$^\circ$, 0$^\circ$), 每层厚度为$t=0.012$, 质量密度$\rho =1.0$. 板中心处受一垂直于板的正弦激励作用$P(t)=-p\sin (200\pi \tau )$, 初始条件为${\dot{{x}}}(0)=0,{x}(0)=0$. 材料属性和正弦激励幅值同样均具有不确定性, 假设其区间范围$I$如表3所示; 同时假设它们在区间上服从正态分布, 均值$\mu $和标准差$\sigma $也如表3所示. 采用4结点矩形薄板单元, 将板划分为相等的16个单元. 为方便, 所有数据均采用无量纲量.

图 5

新窗口打开| 下载原图ZIP| 生成PPT

图 5 四边固支复合材料层合板

Fig. 5 A composite laminate with fully clamped boundary conditions

表 3 材料属性和正弦激励幅值的不确定性

Table 3 Uncertainties of material properties and the amplitude of the harmonic excitation

新窗口打开| 下载CSV

考虑板中心处$z$方向的动力响应. 在同一时间步长$\Delta \tau =1.0\times 10^{-5}$下, 本文所提随机、区间方法和传统随机、区间方法在$\tau =0\sim 0.05$内计算得到的响应范围以及响应标称值曲线如图6所示, 其中图6(a) ~图6(b)分别为两种随机、区间方法得到的响应曲线, 图6(c)为本文所提随机、区间方法计算得到的响应曲线. 本文所提随机、区间方法和传统随机、区间方法计算得到的$\tau=0.025$时刻位移上下界如表4所示.

表 4 本文所提随机、区间方法和传统随机、区间方法计算得到的板中心处$z$方向在$\tau=0.025$时刻的位移上下界

Table 4 The upper and lower bounds on the displacement of the center of the plate in $z$-direction at $\tau =0.025$ obtained by SHPM, IHPM TSM and TIM

新窗口打开| 下载CSV

由图6和表4可知, 本文所提随机、区间方法得到的响应上下界曲线分别与传统随机、区间方法所得响应上下界曲线同样近乎完全重合,再次验证了所提方法的有效性. 同时, 本文所提区间方法得到的响应上界比随机方法得到的响应上界大, 区间方法得到的响应下界比随机方法得到的响应下界小, 也与理论推导和算例6.1现象一致.

图 6

新窗口打开| 下载原图ZIP| 生成PPT

图 6 本文所提方法和传统方法计算得到的板中心处$z$方向的响应曲线

Fig. 6 The response curve of the center of the plate in $z$-direction obtained by SHPM, IHPM, TSM and TIM

下面进一步比较在时间步长$\Delta \tau =2.5\times10^{-4}$下利用哈密顿系统辛算法和直接求解运动微分方程计算得到的响应标称值,曲线分别如图7(a) ~图7(b)所示.

图 7

新窗口打开| 下载原图ZIP| 生成PPT

图 7 时间步长$\Delta \tau =2.5\times 10^{-4}$下利用不同算法计算得到的板中心处$z$方向的响应标称值

Fig. 7 The nominal value of the response of the center of the plate in $z$-direction obtained by different algorithms with time step $\Delta \tau =2.5\times 10^{-4}$

由图7可知, 利用哈密顿系统辛算法求解能够得到满意的响应标称值, 然而在该时间步长下直接求解运动微分方程得到的响应标称值同样在很短时间内发散. 这一现象再次体现了哈密顿系统辛算法由于能够保持体系结构特征而具有很好的稳定性, 凸显了哈密顿系统辛算法在数值模拟中的显著优越性.

7 结论

本文考虑非齐次线性哈密顿系统, 基于摄动理论发展了含小参数扰动的哈密顿系统的参数摄动法, 以此为基础分别将小扰动看作随机、区间扰动,提出了分别针对随机、区间不确定性哈密顿系统的参数摄动法, 突破了传统哈密顿系统的限制, 并由切比雪夫不等式证明了两者响应分析结果的相容性关系, 通过两个综合考虑外载荷和不确定性的结构动力响应数值算例验证得到了如下结论:

(1)在同一较大时间步长下, 由直接求解运动微分方程的方法得到的响应标称值很可能会发散, 而将运动微分方程引入到哈密顿系统后再利用辛算法计算仍然可以得到令人满意的结果, 体现了哈密顿系统辛算法具有很好的稳定性, 保持体系结构特征的优势明显, 在数值仿真模拟方面的显著优越性;

(2)在同一较小时间步长下, 本文所提随机不确定性哈密顿系统的参数摄动法得到的响应上下界分别与传统随机方法所得响应上下界结果近乎完全一致, 差别微小, 本文所提区间方法得到的响应上下界也分别与传统区间方法所得响应上下界极为相近, 验证了所提方法的有效性;

(3)本文所提随机方法在较大时间步长下得到的响应上下界与蒙特卡洛模拟方法在较小时间步长下得到的响应上下界基本吻合,说明本文所提方法在较大时间步长下也能得到具有较高精度的结果, 再次验证了本文所提方法的优越性;

(4)当参数的区间范围由概率统计信息确定, 即区间半径表示为距其均值距离为标准差的整数倍时, 本文所提区间方法计算得到的响应区间范围包含本文所提随机方法计算得到的响应区间范围, 即区间响应上界大于随机响应上界、区间响应下界小于随机响应下界.

参考文献

原文顺序

文献年度倒序

文中引用次数倒序

被引期刊影响因子

[1]

李鸿晶, 梅雨辰, 任永亮 .

一种结构动力时程分析的积分求微方法

力学学报, 2019,51(5):1507-1516