（固体力学专业论文）单自由度非线性随机系统的响应.pdf

摘要 本文提出构造同时受高斯自噪声外激和参激作用下的强非线性单自由度随 机系统精确平稳解的新方法。该方法主要思想是寻找该系统的等价非线性随机 系统的精确平稳解以替代原系统的解。本文得到了求解等价随机系统的方法， 给出了该系统解存在的条件，给出了精确平稳解的形式，并且证明已得到的能 量依赖的精确平稳解只是它的特殊情形，其他几类都是能量非依赖的解，因此 本文得到的精确平稳解是迄今为止最广泛的一类单自由度强非线性随机系统的 精确平稳解。 。 由于只有在严格的限制条件下才能得到高维f p k 方程的精确平稳解，通常 情况下只能寻求其数值解。本文采用有限差分法将f p k 方程在空间域上进行离 散，通过求解微分方程通解得到f p k 方程的瞬态解。针对线性系统和非线性系 统分别给出了几个典型时刻的条件转移概率密度和各分量的响应矩随时间交化 曲线。针对简化的平稳f p k 方程，同样采用差分法结合特征值展开求解其稳态 解，并与差分法结合超松弛迭代法的求解结果进行比较，讨论了差分法的阶次、 网格密度等对数值解精度的影响。该数值解法对于求解f p k 方程有相当高的精 度，因此本文所提数值算法是求f p k 方程近似解的有效路径。 关键词：精确平稳解，非线性随机系统，等效系统， 差分法，瞬态解， f p k 方程 a b s t r a c t i nt h i sp a p e r , an e w p r o c e d u r ei sp r o p o s e dt oo b t a i l lt h ee x a c ts t a t i o n a r ys o l u t i o n so f s o r u ec l a s s e so ft h es i n g i e - d e g r e e - o f - f i e e d o m ( s d o f ) s t r o n g l yn o n l i n e a rs t o c h a s t i c s y s t e m ss u b j e c t e dt op a r a m e t r i ca n d o ra x t e l l l a lc r a u s s i a nw h i t en o i s e s t h em a i n i d e ao ft h ep r o c e x i n r ei st of i n dt h ee x a c ts t a t i o n a r ys o l u t i o no ft h ee q u i v a l e n t n o n l i n e a rs t o c h a s t i cs y s t e m t h ec o n d i t i o n sf o rt h ee x i s t e n c eo f t h ee x a c ts t a t i o n a r y s o l u t i o na r ed e s c r i b e d ，a n dt h ef o r m so ft h ee x a c ts t a t i o n a r ys o l u t i o n sa r eo b t a i n e d i ti ss h o w nt h a tt h ek n o w ne x a c ts t a t i o n a r ys o l u t i o n sw h i c ha r ed e p e n d e n to fe n e r g y a r eo n l yt h es p e c i a lc a s e so fs o l u t i o n sp r e s e n t e di nt h i sp a p e r , w h i l es o m eo t h e r c l a s s e so fs o l u t i o n sa l lp o s s e s st h ep r o p e r t yt h a tt h es o l u t i o n s 锄i n d e p e n d e n to f e n e r g y t h u s t h ee x a c ts t a t i o n a r ys o l u t i o n so b m i n o di nt h i sp a p e ra 托t h o s e b e l o n g i n gt ot h em o s tg e n e r a lc l a s so ft h es d o fs t r o n g l yn o n l i n e a rs t o c h a s t i c s y s t e m ss of a r t h ee x a c ts t a t i o n a r ys o l u t i o n sf o rf o k k e r - p l a n c k - k o l m o g o r o v ( f p k ) e q u a t i o na r e u s u a l l yo b t a i n e du n d e rs t r i c t c o n d i t i o n s t h ef p ke q u a t i o ni s 垭m l l ys o l v e d n u m e r i c a l l y s e v e r a lm e t h o d ss u c ha sp a t hi n t e g r a lm e t h o d f i n i t ed i f f e r e n c em e t h o d e t a 1 h a v eb e e ne m p l o y e dt on u m e r i c a l l ys o l v et h ef p ke q u a t i o n t h ef i n i t e d i f f e r e n c em e t h o di sa p p l i e dt os o l v et h ei s s u eo fs p a t i a ld i s c r e t i z a t i o ni no r d e rt o t r a n s f o r mt h ef p ke q u a t i o nt oas y s t e mo fo r d i n a r yd i f f e r e n t i a le q u a t i o n s t h e nt h e n o n s t a t i o n a r ya n ds t a t i o n a r ys o l u t i o n sa l eo b t a i n e db ys o l v i n gt h es e to fo r d i n a r y d i f f e r e n t i a le q u a t i o n s s e v e r a le x a m p l e sa g i v e nt oi l l u s w a t et h ea p p l i c a t i o no f p r o p o s e dp r o c e d u r e s i ti ss h o w nt h a tt h i sm e t h o di sav e r yg o o dp r o c e d u r ef o r n u m e r i c a ls o l u t i o no ff p k e q u a t i o na n dt h em e t h o dc 锄b ef i a t h e ra p p l i e dt oh i g h e r d i m e n s i o n a lf p k e q u a t i o n k e yw o r d s ：e x a c ts t a t i o n a r ys o l u t i o n , n o n l i n e a rs t o c h a s t i cs y s t e m ， e q u i v a l e n ts y s t e m , f i n i t ed i f f e r e n c em e t h o d ，n o n s t a t i o n a r y s o l u t i o n 独创性声明 本人声明所呈交的学位论文是本人在导师指导下进行的研究工作及取得的 研究成果。据我所知，除了文中特别加以标注和致谢的地方外，论文中不包含其 他人已经发表或撰写过的研究成果。也不包含为获得逝姿盘茔或其他教育机 构的学位或证书而使用过的材料。与我同工作的同志对本研究所做的任何贡献 均已在论文中作了明确的说明并表示谢意。 学位论文作者签名；签字日期：年月日 学位论文版权使用授权书 本学位论文作者完全了解逝姿盘堂有关保留、使用学位论文的规定， 有权保留并向国家有关部门或机构送交论文的复印件和磁盘，允许论文被查阅和 借阅。本人授权滥江态堂可以将学位论文的全部或部分内容编入有关数据库 进行检索，可以采用影印、缩印或扫描等复制手段保存、汇编学位论文。 ( 保密的学位论文在解密后适用本授权书) 学位论文作者签名： 签字日期：年月日 学位论文作者毕业后去向： 工作单位： 通讯地址： 一名：墨锄 签字日期：叼年月护日 电话： 邮编： 浙江大学硕士学位论文 第一章绪论 在自然及社会领域中存在着大量的随机扰动，如地震、风、浪波、喷气噪声 等等，乃至股市的波动及社会的动荡等等都是非常常见的随机扰动。关于随机动 力学的研究最早可上溯至上世纪初爱因斯坦关于布朗运动的研究，由于上世纪五 六十年代航空航天等方面的需求，诞生了随机振动这一学科。开始，大多数随机 振动理论研究都基于确定性时不变线性动态系统模型，很快就建立了较为完善的 线性随机振动理论。此后，许多研究者的注意力转向了非线性随机振动问题，发 展了多种预测非线性随机响应的方法，部分地揭示了随机振动中的非线性现象。 随机系统的稳定性与参激随机振动理论也有较大地发展。同时还发展了结构对宽 带随机激励的响应的渐近分析方法，揭示了结构均方响应的空间分布的一些渐近 规律。在随机振动的线性与非线性动态系统的可靠性方面也提出了多种近似方 法。随着随机振动理论与技术的发展，其应用也越来越广泛。开始阶段主要应用 于运载工具( 飞行器、汽车、船舰等) ，后来扩展到各种复杂的结构( 高层建筑、 离岸结构等) 。今天，随机振动已发展为一门内容十分丰富的学科【l 捌。 对非线性随机振动理论的研究经过几十年的研究，虽然已发展了许多方法， 但总的来说尚很不成熟。对于一扩散的马尔可夫过程，其响应的转移概率密度可 用福克普朗克一柯尔莫哥洛夫( f p k ) 方程描述， 设r l 维时不变随机动力系统，由以下随机微分方程描述 d x l = q ( x ) d t + ( x ) d b ( t ) i = l ，月；甜= 1 ，z 其中q 为漂移系数，= 吒o 为扩散系数，e o ) 为独立并具有单位强度维纳过 程。 与肺随机微分方程( 1 - 1 ) 相应扩散过程的转移概率密度满足如下的f p k 方程 垒：一旦业+ 三业 d t 瓠| 2 瓠i 瓠i d ep = p ( x ，t l , o ，t 。) 为条件转移概率密度，并满足如下初始条件 p = p ( x , t l x o ，f 。) i ，嘞= 6 0 x o ) 1 ( 1 - 2 ) ( 1 - 3 ) 浙江大学硕士学位论文 转移概率密度还需满足如下归一化条件 i p ( x ，r k ，t o ) d r = l ( 卜4 ) 求该f p k 方程的解就得到该系统响应的统计分布。已得到f p k 方程精确瞬态 解的主要是受高斯白噪声外激的时不变线性系统，其解最早由w a n g 与 u h l e n b e c k 【3 】用傅立叶变换技术得到，它也可用特征函数展开法求得。非线性系 统f p k 方程的精确解非常难得到，迄今已得到精确解的只是一些一阶非线性系 统及多自由度常系数线性系统【5 - 7 1 。 目前已碍到的较多的解是精确平稳解，所谓精确平稳解是满足不含时间导数 项的简化f p k 方程的解。当时间足够长时，条件转移概率密度趋于稳定且与初 始条件无关，其稳态概率密度满足的简化f p k 方程为 三翌鱼盟一旦业：o ( 1 - 5 ) 2 融t 。x j e x i i ，j = l ，2 ，n 其中p ( 而，矗) 为稳态概率密度且满足归一化条件 它鼍 i ip ( 五，) 幽呶= l ( 卜6 ) 几十年来，科学工作者一直对非线性随机系统的响应感兴趣。二阶非线性随 机系统只能由平稳f p k 方程得到精确平稳解。首先a n d r o n o v 等【8 】与k r a m e r s l 9 独立地得到在g a u s s 白噪声外激下非线性弹簧系统的精确平稳解。c a u g h e y 与 p a y n e 1 0 】及c a u g h e y 与m a 1 1 , 1 2 将该解推广于g a u s s 白噪声外激下某些类型的非线 性阻尼系统。d i m e n t b e r g 1 3 】首先得到同时受g a u s s 白噪声外激和参激的特殊单 自由度系统的精确平稳解。y o n g 与l i n 1 卅用详细平衡方法得到了一大类同时受 g a u s s 白噪声外激和参激的非线性随机系统。l i n 与c a i 1 5 1 提出了广义平稳势法 得到非线性随机系统的精确平稳解。z h u ( 朱) 1 6 1 将c a u g h e y 与m a 1 1 , 1 2 及 d i m e n t b e r g ”1 之解推广于更为一般的几类非线性随机系统。 为研究多自由度非线性随机系统的精确平稳解，f u l l e r 1 7 1 将在g a u s s 白噪 声外激下的非线性系统表示成随机激励的耗散的h a m i l t o n 系统的形式，并将相 应f p k 方程的平稳解表示成h a m i l t o n 函数的泛函。8 0 年代末与9 0 年代初， s o i z e 1 砌h 与z h u 等【2 1 2 2 独立地得到了g a u s s 白噪声外激与参激下相当一般的非 2 浙江大学硕士学位论文 线性耗散h a m i l t o n 系统的f p k 方程精确平稳解。所有这些精确平稳解都具有能 量等分的性质，即所有这些平稳解都是能量的泛函，能量等分解意味着一个多自 由度非线性随机动力学系统的各自由度之间的能量比值是固定的，改变随机激励 的大小，只能改变系统总能量的分布，而不能改变各自由度间的能量分配。然而 g a u s s 白噪声外激下的时不变线性系统的平稳概率密度是g a u s s 的，改变随机激 励与阻尼的大小，可同时改变系统总能量的分布与各自由度间能量的分配，即它 不是能量等分的，因此在线性与非线性多自由度随机动力学系统的精确平稳解之 间存在之不一致性。为了克服这一不一致性，近年来，z h u 、h u a n g 等 2 3 - 2 5 1 通过 将g a u s s 白噪声激励的多自由度非线性系统由l a n g r a n g e 体系转化为h a m i l t o n 体系，利用h a m i l t o n 系统的可积性与共振性，得到了完全不可积、完全可积非 共振、完全可积共振、部分可积非共振、部分可积共振五类随机激励的耗散的 h a m i l t o n 系统的精确平稳解，其中第一类是能量等分的，后四类是能量非等分的， 从而打破了一直以来只有非线性随机系统能量等分精确平稳解的局面。文1 2 6 1 则 得到了最一般的p o i s s o n 白噪声激励下的精确平稳解。 能求精确瞬态解的f p k 方程相对来说是少数，大多数工程实际问题中，往往 只能求近似解析解或数值解，对高阶f p k 方程的数值解法是随机振动理论走向实 际应用必须要解决的重大课题。现已发展了一些数值求解方法，如蒙特卡罗模拟 法、随机步行法、路径积分法、有限元法、有限差分法等，下面分别作简要介绍。 1 1 蒙特卡罗模拟法”刀 蒙特卡罗模拟法又称数字模拟或随机模拟，该方法的基本思想是概率统计， 例如可用抛硬币的方法得到随机样本。在随机振动中，蒙特卡罗模拟法利用 计算机产生的随机样本来近似求解系统响应的统计特性，近似程度随随机样本数 的增加而提高，蒙特卡罗模拟法有如下三个基本步骤 i 、由于计算机真正能产生的随机数只有0 、1 两个数，即高电平和低电平， 要根据具体的系统激励过程把计算机原始的随机数转化为所需要的激励样本，由 算法转化得到的随机激励样本要与系统激励的概率统计特性尽量保持一致，算法 要尽量简练，这样才能提高蒙特卡罗模拟法的计算速度。 浙江大学硕士学位论文 2 、将随机激励样本代入随机运动方程求解其响应样本。 3 、从大量响应样本求解其统计特性，如概率密度、谱密度等。 蒙特卡罗模拟法是在随机振动中应用最广的数值计算方法，它是从原方程 ( 1 - 1 ) 出发进行统计分析的，适用于高维的线性与非线性随机系统，所以在随 机振动中是应用最广的数值计算方法，而且不受计算机存储误差的干扰，其精度 也比较高。蒙特卡罗模拟法的另一个优点是所产生的各个随机样本及求解得到的 响应样本相互之间没有联系，非常适合于用并行算法进行计算。随着计算机性能 的不断提高，用蒙特卡罗模拟法结合并行算法对高维非线性随机系统的响应计算 范围将不断扩大。当然该方法也有缺点，其主要缺点是提高精度困难，蒙特卡罗 模拟法所得到的结果精度受激励样本的约束而带有随机性，且结果统计的误差随 样本数( m ) 的增大而减小，其误差a s 1 d m ，因此要增多一位有效数字样 数本要增大l o o 倍。为提高样本的统计精度，特别是小概率事件的统计精度，又 提出不少的改进统计方法，如重要性样本方法等。 1 2 随机步行法啪3 卵 t o l a n d l l 9 1 ，r o b e r t s l 2 0 1 等人曾经用随机步行求解过f p k 方程。如图1 1 所 示，随机步行是指一个人在离散时刻f ，= ，& ( ，= 0 ，1 ，2 ) 上沿y 轴随机地向 左或右走一步，人所在的位置只能是离散值y k = k a y ( k = 0 ，士1 ，2 ) ，这是必 然的任何时刻人所在的所有位置不可能离开雎。假设时间t y = 皿时人在 y k = k a y 上，当f ，+ l = u + 1 ) a t 时刻，人向y k “= + 1 ) 缈走的概率设为，折人 向以。= ( k - d a y 走的概率为卜既。随机步行是一种概率扩散过程，其一二阶 增量矩阵分别是 ， a a ：警：篡a t z ：鼍。骈：。a y ) ： 。 c ，们 2 饥，0 ，) = 岛( 劬2 + ( 缈) 2 = ( 2 为了与一维f p k 方程的概率进化过程相匹配，要使 4 浙江大学硕士学位论文 魄) 砂= a l ( y k ，t _ ，) a t ( 缈) 2 = 6 1 l ( 儿，t , ) a t ( 1 - 8 ) 其中巩+ - - 1 ，o l ，岛。为漂移与扩散系数，可由方程( 卜1 ) 给出。在给定缈后， 由漂移系数与扩散系数q 饥，) ，6 1 。魄，) 得到岛，和r 。 概概率钆 夕 y k 4y ky h l 图1 1 随机步行法的示意图 用随机步行法求解f p k 方程的学者不是很多，其原因是该算法不是无条件稳 定的，当所设的时间步长大于它的条件稳定的最大步长时将使该随机步行法变得 不稳定。因此，用此方法求解时要用较小得时间步长，且该法很难推广于高维 f p k 方程得求解。 1 3 路径积分法“州1 1 路径积分法( p a t h - - i n t e g r a lm e t h o d ) 是一种可用于求解较高维f p k 方程 的数值方法，它的基本思想是将f p k 方程之解表示为如下形式的泛函积分 丘 p ( y ，r ) = i d u ( y ) e x p - i 考 穸( r ，) ，y ( f ，) 弦7 ) p ( y o ，) ( 1 9 ) y o 如 式中d u ( v ) 是一个积分测度， 称为o n s a g e r - m a e h l u p 泛函。在空间与时间上离 散化，以路径和代替路径积分，则 5 浙江大学硕士学位论文 1一nn-i p ( y ，f ) - j 罂骡m 呐e x p - r 蔡( y k + 1 , y k f ) ) p ( y )呻 ，2 u ”+ 吨 ( 卜1 0 ) = ，。j i mn 协，g ( y “，y t 彳) p c v o , t ) n n r - 1 - t o k j 式中g k + l ，y i f ) = u ke x p - - o , k = 0 , 1 , 2 , 3 ) ，由方程( 2 - 1 3 ) 得 ，“，而) = 一 兰+ 眈2 ( a o + a l x l 2 + a 2 2 2 + # ) 斗? 2( 2 硝) 丝古蔓塾五蔓z 竖! 二兰旦f 竺2 兰堡3 i 兰2 2 ( 口o + 口l # + 口2 + a 3 群) 该系统的稳态概率密度为： p ( 而，毛) ；c ，( a 。+ 口。彳+ 口：+ a ，# 弦( d “埔嗨考+ 畸彳弓每+ 扣奇j ( 2 - 4 5 ) 其中c 为归一化常数。 例6 ； 浙江大学硕士学位论文 令c ( 而，艺) = 一芹”2 ”衙引， 女2 0 ，七= 1 ，2 ，3 ) ，由方程( 2 1 3 ) 得 胍班丢嚣+ 町慨札扣膏珈+ 脚：鼎) + 去b 州硝+ 靠h 1 觚妒枷撕删 该系统的稳态概率密度为： 肌。、一r 十”加蒯小( 椭去一砷啪 p ( 而，恐) = c p o 吒”何j 2 3 本章结论 ( 2 - 4 6 ) ( 2 - 4 7 ) 本章提出了求解受高斯白噪声外激和参激的单自由度强非线性随机系统的 精确平稳解的方法。它本质上是通过求解等效随机系统f p k 方程的解代替原系 统的解，我们得到了寻找该等效随机系统的条件及求解方法，并由该方法得到了 几类强非线性随机系统的精确平稳解，得到的精确平稳解的形式不仅包含了能量 依赖的解还有能量非依赖的解。对于能量依赖的解包括了已有最一般的强非线性 系统精确平稳解，而其他的能量非依赖系统的解是全新的解，该方法还可能推广 于求多自由度强非线性随机系统的精确平稳解。 2 0 浙江大学硕士学位论文 第三章用差分法求f p k 方程的数值解 3 1 前言 如第一章所述，f p k 的精确解很难得到，为了求其瞬态解， s f ：w o m 蚀、1 c z ，l a b e r ( 3 m a n ，b e s p e n c e r 等在文【6 5 】中采用有限单元 和有限差分法，同时对时间和空间进行离散，得到线性代数方程组，从而可以求 解每一时刻的概率密度分布。 即便是精确平稳解也是极其有限的且只有在非常严格的条件下才能得到。虽 然在第二章我们对某几类强非线性系统，可以用等效非线性系统法及随机平均法 等求得其平稳解以扩大其应用范围，但也需要在一定的条件下才能得到。在通常 情形下需数值求解该f p k 方程，国内学者张森文【6 6 j 用小波法求解了f p k 方程的 平稳解，林家浩【6 7 】等提出用虚拟激励法求解平稳与非平稳随机振动系统的解。 这些方法在研究某些特定系统的近似解方面取得了成功，但至今缺少一种高效简 洁通用的直接计算高维平稳f p k 方程的有效算法，高阶f p k 方程的数值解法是随 机振动理论走向实际应用必须要解决的重大课题。 差分法是求解高维偏微分方程的有效方法之一，而高维稳态f p k 方程的解往 往是某局部区域的一种分布，当边界条件取得足够大时，边界效应影响很小，这 刚好避免了差分法对边界条件比较敏感的缺点且不存在不规则网格，因此用差分 法求解高维f p k 方程的稳态解是很有前途的。文【缱】用差分法直接求解t - - 维f p k 方程的稳态解，研究了不同的差分格式得到近似解的精度。而直接用数值方法求 解高维f p k 方程的最大困难是随着求解维数的提高，代数方程组的未知量以几何 级数增长，例如求解4 维的稳态f p k 方程，求解区域按4 0 x 4 0 4 0 4 0 均匀网 格离散化，未知量个数为4 0 4 ，用二阶差分格式1 6 9 需要用超级计算机及并行处理 方法求解4 维f p k 方程的近似解。 本章以差分法为基础来求解f p k 方程的瞬态解和精确平稳解，通过单自由度 线性与非线性系统的几个算例说明所采用的数值解法的可行性。 浙江大学硕士学位论文 3 2 用差分法求f p k 方程瞬态解 f p k 方程( 1 - 2 ) 只有在特定条件下才存在精确瞬态解，一般情形下只能数 值求解，现用差分法进行求解。此处采用二阶差分格式。 利用方程( 1 - 2 1 a ) 的差分格式将f p k 方程( 1 - 2 ) 在空间上进行离散，可 得线性微分方程组， 【k 】 p ) = 磊d p ， ( 3 1 ) 其中系数矩阵 k 】与不同差分格式及不同的系统特征有关，p = ( p i ，) 7 为空 间离散化后各节点的转移概率密度组成的矢量，且各点的转移概率密度还需要满 足归一化条件。 求解系数矩阵 k 】的特征值和特征向量，由常微分方程的知识可知，方程组 ( 3 1 ) 的通解为， p ( f ) = c i 巧 ( 3 2 ) 其中，q 为任意复常数，其值由初始条件确定； 为特征值( 一般为复数) ，k 是 对应于特征值凡的非零解( 特征向量) ，由( 3 2 ) 式可得到任意时刻的概率密度 分布。 算例1 考虑受高斯白噪声外激的线性系统 舅+ 卢量+ 嘞2 x = w ( t ) ( 3 - 3 ) 其中w ( t ) 为强度2 d 的高斯白噪声。设y = 章，与( 3 3 ) 式相应的f p k 方程为 曼丝丝一a ( y p ) 一盟! 二墅二堡生田：塑 钞2 钞 西 ( 3 - 4 ) 其中p = p ( 墨j ，) 为位移与速度联合转移概率密度；边界条件为，一o ， h ，i y l 寸m ；初始条件为烈而弘o ) = 6 ( y 一0 7 5 ) 。 当卢= o 1 ，d = 0 0 1 ，= 1 时，由前所述求解瞬态解的方法可得f p k 方程( 3 4 ) 浙江大学硕士学位论文 f p k 方程存在精确瞬态解。精确解的表达式为渺7 1 】 p ( x ，f ) - - - ( 2 ，r ) ”| r ( r ) i 】1 ”e x p - 吉 x - u ( t ) 7 r 。( f ) 【z p ( r ) 】 ( 3 5 ) 其中，x = 卜y 】。，p ( f ) 和r ( r ) 分别表示矢量x ( f ) 的均值和协方差矩阵。对于 二维系统，均值矢量和协方差矩阵可准确地表示为时间的函数， m = 去 纂二鬈嚣：飘 净e ， r ( f ) = 而2 d 其中，口。= 譬+ 业2f ，口：= 壁2 一塑奎2 生j ，x o 为初始时刻x 、y 的 值。 由于初始条件取6 函数非常困难且可能引起数值不稳定性，为了避免这种不 稳定性带来的求解误差，我们将t o = l o s 时的转移概率密度作为新的时间起点 f 0 = 0 ，以此时的概率分布作为初始条件的概率分布。采用2 阶差分格式进行离散， 4 0 4 0 网格密度。 图3 卜3 4 分别给出了几个典型时刻的转移概率密度分布，其中( a ) 为由 式( 3 5 ) 得到的精确瞬态解，( b ) 为用数值方法得到的解。 图3 5 3 9 分别为位移均值、速度均值、位移方差、速度方差、位移与速度 相关矩随时间t 的变化。由图3 1 - 3 9 知，对于单自由度线性随机系统，数值解 与精确瞬态解相比有较高的精度。 毒 以 幽塑箐坐 堡 竺2萃 兰敏矿 浙江大学硕士学位论文 1 5- 1 5 ( a ) 精确解 - 1 51 5 ( b ) 数值解 图3 1t - - o ( s e e ) 时刻的转移概率密度 浙江大学硕士学位论文 1 5 - 1 5 ( a ) 精确解 1 5- 1 5 ( ”数值解 图3 2 仁1 0 ( s e c ) 时刻的转移概率密度 浙江大学硕士学位论文 - 1 5- 1 5 ( a ) 精确解 - 1 5- 1 5 ( b ) 数值解 图3 3 仁2 0 ( s e e ) 时刻的转移概率密度 x ，y ) 堑垩查兰堕主兰堡丝苎 - 1 51 5 ( a ) 精确解 1 5 5 ( b ) 数值解 图3 4t = - 3 0 ( s e c ) 时刻的转移概率密度 浙江大学硕士学位论文 o 5 0 4 o 3 o 2 g o _ 1 山0 o 1 o 2 - o 3 - 0 4 o123 4567891 0 t ( s e c ) 图3 5 位移均值随时间的变化曲线 ( 数值解；一精确解) o12345e789 o t ( s e c ) 图3 6 速度均值随时间的变化曲线 ( 数值解；一精确解) 5 4 3 2 1 o f 2 3 4 o o o o n 也 也 也 也 苫山 浙江大学硕士学位论文 o1234567891 0 t ( s e c ) 图3 7 位移方差随时间的变化曲线 ( 数值解；一精确解) o 123 4 5b7891 0 t ( s e c ) 图3 8 速度方差随时问的变化曲线 ( 数值解；一精确解) 舶 埘 忽 妣 m = 罢 m m m m o o o o o o o o o f x l 哪 m 毖 侣 = 2 住 啦 眦 m 玑 f a 山 浙江大学硕士学位论文 图3 9 位移和速度相关矩随时间的变化曲线 ( 数值解；一精确解) 算例2 考虑受高斯白噪声激励的杜芬振子，其运动方程为 2 + p i 一2 x + s x 3 = w ( t ) ( 3 - 8 ) 其中w ( t ) 为强度2 d 的高斯白噪声。设y = 圣，与( 3 8 ) 式相应的f p k 方程为 坳一塑! 一a ( - s y + a , o ：x - s x ) p ：至( 3 - 9 ) 勿2 o x 砂 a 其中p = p k 舅f l 而，t o ) 为位移与速度联合转移概率密度，边界条件为p o ， i 卅，酬啼m ；设初始条件为正态分布 p ( x o , y o , o ) = 2 - - - e x p t 掣+ 学， 对不同的初始分布，其数值瞬态解是不同的，f i g z ，= o ，几= 0 ，- - o y 2 = o 0 5 时为初始条件1 ，取以= o ，以= o 2 5 ，一- - u ，2 = o 0 5 时为初始条件2 均采用 2 阶差分格式进行离散，4 0 x 4 0 网格密度。 浙江大学硕士学位论文 图3 1 0 3 1 9 给出了几个典型时刻的概率密度分布，其中( a ) 、( b ) 分别为初始条 件1 和2 时的概率密度分布。通过对比，可以看出初始条件取为1 时方程很快趋 于稳定，其收敛速度大于取初始条件2 时的情况。图3 2 0 一3 2 4 给出了位移均值、 速度均值、位移方差、速度方差、位移和速度相关矩随时间的变化 1 5- 1 5 ( a ) 初始条件1 1 5- 1 5 ( b )初始条件2 图3 1 0 t - - o ( s e e ) 时的联合概率分布 3 l y 激江大学硕士学位论文 p ( x 。 - 1 5- 1 5 ( a )初始条件l - 1 5- 1 5 ( ”初始条件2 图3 1 1 卢1 ( s e c ) 时的联合概率分布 y 搬江大学硕士学位论文 - 1 51 5 ( a )初始条件1 - 1 5- 1 5 ( b )初始条件2 图3 1 2t = 2 ( s e e ) 时的联合概率分布 浙江大学硕士学位论文 p ( x ， - 1 5- 1 5 ( a ) 初始条件1 x ，y ) - 1 5- 1 5 ( ”初始条件2 图3 1 3t = 3 ( s e c ) 时的联合概率分布 浙江大学硕士学位论文 1 51 5 ( a )初始条件1 - 1 5- 1 5 ( b ) 初始条件2 图3 1 4t - - 4 ( s e e ) 时的联合概率分布 浙江大学硕士学位论文 - 1 5- 1 5 ( a )初始条件1 x ，y ) - 1 51 5 ( ”初始条件2 图3 1 5 t = 5 ( s e e ) 时的联合概率分布 浙江大学硕士学位论文 1 5- 1 5 ( a )初始条件1 - 1 5- 1 5 ( ”初始条件2 图3 1 6t = l o ( s e e ) 时的联合概率分布 3 7 塑垩查兰堡主兰垡丝塞一 - 1 5- 1 5 ( a ) 初始条件1 - 1 5- 1 6 ( b )初始条件2 图3 1 7t - t - t 兑o ( s e c ) 时的联合概率分布 浙江大学硕士学位论文 5- 1 5 ( a )初始条件1 - 1 5- 1 5 ( b )初始条件2 图3 1 8t = 3 0 ( s e c ) 时的联合概率分布 x ，y l 浙江大学硕士学位论文 - 1 51 5 ( a )初始条件1 - 1 5- 1 5 ( b ) 初始条件2 图3 1 9t - - 4 0 ( s e e ) 时的联合概率分布 4 0 x ，y ) 浙江大学硕士学位论文 o 0 3 o 0 2 o 0 1 o o 0 1 - 冒- - o 0 2 - o 0 3 - o 0 4 - o 0 5 - o 0 6 o 0 7 o 3 o 2 5 o 2 一o 1 5 u o 1 o 0 5 o o 0 5 - l 倒瑁系佣：1 f 厂 ， 厂一 j l ， ， 一 初穿日缘午z | ； v o51 01 52 02 53 0 t ( s e c ) 图3 2 0 位移均值随时间的变化曲线 。 7 f 。l r 。 ，初始嗣件2 ， ， - a ：r t 馅z【埘! ， - _ - - - o51 01 52 02 5 3 0 t ( s e c ) 图3 2 1 速度均值随时间的变化曲线 4 l 浙江大学硕士学位论文 o 1 1 o 1 o 0 9 o o 0 5 o 2 5 o 1 5 o 1 l 一 。 ，一 | ， 弋 l f 7 | 初j眙条干f 1 _ 。仞始爰酷件2 i o12345678 9 1 0 t ( s e c ) 图3 2 2 位移方差随时间的变化曲线 厂 |厂、1宅 】始条件1 y 一一 。每 - - - - - i 太= 多 j 。 入 t芴始条件2 i7 o1234567891 0 t ( s e c ) 一 图3 2 3 速度方差随时间的变化曲线 o 【z x l 山 譬山 浙江大学硕士学位论文 o 1 o 0 8 o 0 0 4 蛊 i l l o 0 2 o o 0 2 - 0 0 4 i 初始j荼件1 | 一l | 1f ，、 i l | |j| 电9 石迤 5 形 、 ：；_ ， ； t 乡 刚厢幂什 o1 2 34567891 0 t s e c ) 图3 2 4 位移和速度相关矩随时间的变化曲线 3 3 用差分法求f p k 方程稳态解 当时间足够长时，条件转移概率密度趋于稳定且与初始条件无关，其稳态概 率密度满足简化f p k 方程( 卜5 ) 及归一化条件( 卜6 ) 。 采用前面求解微分方程通解的方法，由线性代数的知识可知，零特征值所对 应的特征向量即是对应f p k 方程的稳态解。 可以采用另一种方法来求解简化f p k 方程( 卜5 ) 的稳态解【弼。由于离散后 的代数方程中，基本带宽内的非零元是很少的，例如应用二阶差分格式，对于二 维的随机系统，f p k 方程经离散化后一个对角元相应的基本带宽内仅有5 个非零 元，该稀疏的代数方程组非常适合用迭代法求解，这可使方程的求解规模大大提 高。因此，用差分法结合迭代法求解高维f p k 方程的稳态解是很有前途的。该方 法不仅计算效率高且有很高的计算精度，是求解高维平稳f p k 方程的有效方法。 同样利用差分格式将平稳f p k 方程( 卜5 ) 离散化，可得线性代数方程组为 k i p l = 0 ( 3 - 1 1 ) 浙江大学硕士学位论文 由于系数矩阵 k 】一般是非正定的矩阵，p = ( a ，岛) 为离散后各节点的 稳态概率密度，且要保证p 满足归一化条件( 卜6 ) ，我们用超松弛迭代法求解方 程( 3 1 1 ) ，即用如下迭代公式6 2 】 ”1 ) = 帕+ 口( 一i - i 巧。“，- 壹巧_ ) 乞c。一12)lm - t ”9 = 帕+ 口l - 巧0 “_ 巧_ i 肛 ( 3 一 其中上标n 表示第n 步迭代的结果，a 为超松弛迭代因子，为系数矩阵【k 】的 元素，为保证式( 3 1 2 ) 的收敛性，口取值范围一般为( o ，2 ) ，其最佳值取决于 矩阵k 的性质，可按下式取值嗽l a 掣2 意萧 1 3 其中p ( k ) 为矩阵k 的条件数。 实际迭代时，初值可取解空间均匀分布，为保证p 的非负性及迭代的收敛性， 将迭代过程中将可能产生的负n 赋为零，且每经1 0 步迭代进行方程( 卜6 ) 的归 一化处理。 由于差分法得到的以节点自由度为基本变量的代数方程组，其带宽内仅有非 常少的非零元，用传统的直接求解法这些带宽内零元素都需要存储单元，当求解 的系统维数大于4 时，数据存储量非常大，例如对4 维随机系统划分为4 0 x 4 0 x 4 0 x 4 0 的网格，则该系统共有4 0 4 自由度，对于二阶差分格式的半带宽约为 4 0 3 ，但在4 0 3 半带宽内实际的非零元素仅为9 个，因此用迭代法求解该方程组对 存储量盼要求大大降低，将使解题规模得到最大限度的扩大。 本节以单自由度非线性随机系统算例对两种以差分格式为基础的为迭代法 及特征值展开法求得的数值解的精度做了比较，讨论了不同网格密度、不同差分 阶次对迭代解精度的影响。 算例 考虑运动方程( 3 8 ) 式的受高斯白噪声激励的杜芬振子，设y = y c ，与其相 应的平稳f p k 方程为 浙江大学硕士学位论文 0 2 ( d 。p ) 一亟型一出生：罂：三二型臼：o ( 3 - 1 4 ) 砂 良 钞 其中p = p ( x ，y ) 为位移与速度联合稳态概率密度，平稳f p k 方程( 3 1 4 ) 式的精 确平稳解为： p c 咖c 唧m 手一争+ 制 洚 其中c 为归一化常数。 当口= o 1 ，d = o o l ，m o = 1 ，艿= 4 时，由式( 3 1 5 ) 得到的位移与速度联合稳 态概率密度精确平稳解及由方程( 3 1 4 ) 用方法1 ：即4 阶差分法及特征值展开 后取零特征值对应的特征向量得到的数值解；用方法2 ：即4 阶差分法及超松弛 迭代得到的数值解分别见图3 2 5 ( a ) ( b ) ( c ) 。本算例数值解用8 0 x 8 0 网格密 度，方法二中超松弛迭代因子口= o 0 6 ，经1 2 0 0 0 步迭代得到。由图3 2 5 知， 方法i 与方法2 都可得到相当精度的精确平稳解。 图3 2 6 ( a ) ( b ) 分别给出了用不同阶的差分格式及网格密度得到的边际概 率密度p ( x ) 与精确解的差的绝对值，采用第一种数值解法。由图3 2 6 知，差分 格式的阶次越高，网格越细，数值解精度越高。图3 2 7 给出了4 阶差分格式时 分别采用方法1 和方法2 得到的边际概率密度p ( z ) 与精确解的差的绝对值。图 3 2 8 给出了4 0 * 4 0 网格密度时分别采用方法1 和方法2 得到的边际概率密度 p ( x ) 与精确解的差的绝对值。可见方法1 与方法2 都能得到很高精度的数值解， 且方法1 比方法2 误差更小。 。 由于方程( 3 - 8 ) 相对应的扩散矩阵是奇异的，因此式( 3 1 4 ) 由差分法得 到的线性代数方程组( 3 1 1 ) 的矩阵k 并不是对角占优的，条件数p ( k ) 较大，因 此收敛的超松弛因子a 很小，本例中口的取值范围为0 0 2 - - 0 0 8 ，在可能的超松 弛因子口的取值范围内，a 取值越大，收敛速度越快。 图3 2 9 给出了差分法与直接由方程( 3 - 8 ) 用蒙特卡罗模拟及用路径积分 法【鹋】得到的数值解与精确解误差的比较，其中蒙特卡罗模拟时用1 0 0 0 个样本， 每个样本用有效的1 0 0 0 0 点进行统计，路径积分法及4