绘制Duffing振子地分叉图地程序
Duffing_vanderPol振子随机分岔的全局分析
Duffing -van der Pol 振子随机分岔的全局分析*徐 伟1)贺 群1)2)戎海武3)方 同4)1)(西北工业大学应用数学系,西安 710072)2)(武警工程学院应用数学系,西安 710086)3)(佛山科技学院应用数学系,佛山 528000)4)(西北工业大学振动工程研究所,西安 710072)(2002年9月20日收到;2002年11月18日收到修改稿)应用广义胞映射方法研究了参激和外激共同作用的Du ffing -van der Pol 振子的随机分岔.以系统参数通过某一临界值时,如果系统的随机吸引子或随机鞍的形态发生突然变化,则认为系统发生随机分岔为定义,分析了参激强度和外激强度的变化对于随机分岔的影响.揭示了随机分岔的发生主要是由于系统的随机吸引子与系统的随机鞍碰撞产生的.分析表明,广义胞映射方法是分析随机分岔的有力工具,这种全局分析方法可以清晰地给出随机分岔的发生和发展.关键词:随机分岔,全局分析,广义胞映射方法,随机吸引子,随机鞍PACC :0547,0545*国家自然科学基金(批准号:10072049)资助的课题.1 引言由于随机噪声干扰在实际系统中总是存在的,考虑随机噪声的影响,特别是对分岔现象的影响成为人们关心的主题之一,随机分岔的研究尚处于起步阶段,如何给出随机分岔的合理定义是研究的核心问题之一.目前的定义主要分为两大类,一类是基于系统的稳态概率密度形状,随参数的变化突然发生变化给出的定义,例如从单峰突然变为双峰,这类定义称为P 分岔[1];另一类是基于系统的最大Lya -punov 指数随参数的符号变化给出的定义,这类定义称为D 分岔[1].研究表明,这两类定义给出的结果并不完全相同.Baxendale [2]给出了一个例子,当Lya -punov 指数符号发生变化时,系统的稳态概率密度形状不依赖分岔参数;另一方面,Crauel 和Flandoli [3]展示了一个相反的例子,当系统的稳态概率密度由单峰变为双峰时,系统的Lyapunov 指数符号不发生变化.究竟如何选择合适的定义展示随机系统的拓扑性质的变化,如何确切的描述随机分岔的发生,这是目前研究的焦点.此外,Meunuer [4],Arnold [1]研究表明,对于一类系统,稳态概率密度的形状变化与随机干扰的变化联系并不密切.这种基于稳态概率密度形状变化的定义,有时很难说明系统的拓扑性质,在噪声干扰下是否发生了真正意义上的变化.而对于Lyapunov 指数的精确计算仍是一个困难的问题,目前虽然有一些算法,例如Wolf [5],但精度并不理想,这直接影响到最大Lyapunov 指数是否为零的判断.这两类定义,均要求系统具有遍历性,而许多系统不满足这一要求,仅能进行局部分析.目前的文献,主要以平均方法进行讨论,通过平均方程的研究,得出原系统的有关结论,但对于许多系统,特别是一类非自治系统,平均方程与原系统的等价性无法保证,这时得出的关于平均方程的有关结论多少能反映原系统的性态,让人产生怀疑.此外,一般的数值计算,采用的是点映射方法,总要进行人为的有限步截断,这在许多讨论中,使得我们对于何时达到稳态产生异议,多少步可以认为是达到了最终的稳态,这对于稳态概率密度的讨论,对于随机跃迁的讨论至关重要.综上所述,对于随机分岔的讨论,仍然存在许多困难和问题,既要表征系统轨线的拓扑性质,又能体现系统的随机特性,确实很困难.我们认为随机分岔定义的关键是如何表征随机系统的拓扑性质,找到随机系统的合适的相对不变量.本文提出以系统的随机吸引子(包括随机鞍)的形态(包括大小、尺寸、周期等)的突然变化表征随机系统的拓扑性质变化,第52卷第6期2003年6月1000-3290 2003 52(06) 1365-07物 理 学 报AC TA PHYSIC A SINICAVol.52,No.6,June,20032003Chin.Phys.Soc.用来描述随机系统的分岔.Duffing-van der Pol振子是一个重要的模型,许多现象经过简化和模型化后可以用这一模型描述.对于确定情形分岔的研究,目前已有不少成熟的结果.Holmes和Rand[6],Guckenheimer和Holmes[7]针对这一模型,讨论了Pitchfork分岔和Hopf分岔,给出了全局分析和局部分析.对于噪声激励的Duffing-van der Pol振子的研究,目前已有一些成果,主要是利用随机平均法,随机规范形方法,在P分岔和D 分岔的定义下,讨论随机分岔的局部分析,这些结果参见Na machchivaya[8]Arnold[1]等,Schenk-Hoppe[9]等人的工作.2 广义胞映射图论方法胞映射的基本思想是把动力系统的状态空间离散化为大量状态胞的集合,用胞映射来描述原系统的动力学行为.Hsu提出两类胞映射方法,一种称为简单胞映射(SC M),一种称为广义胞映射(GC M).广义胞映射与简单胞映射相比,每个胞可以具有多个象胞.广义胞映射方法将原动力系统对应于一个胞映射动力系统.对于一个胞映射动力系统的分析方法可分为两个发展阶段,在早期的分析过程中,Hsu 将胞映射动力系统等价于一个有限时齐的马尔可夫链.采用马尔可夫链的分析方法对相应的胞映射动力系统进行分析,这种分析方法可以计算出相应动力系统的吸引子和吸引域,边界集合,确定吸引子的极限概率分布,以及奇怪吸引子的空间概率特征.所以广义胞映射是研究随机动力系统自然而有力的工具.Guder和Kreuzer[10]证明在一定条件下,胞映射的吸引子、吸引域和边界收敛于原动力系统的吸引子、吸引域和边界.从1995年开始,Hsu[11]在以前的分析方法的基础上,借助于图论方法,对胞映射动力系统的分析方法进一步发展.突出了对系统的拓扑特性的分析.对于胞映射系统,将每个胞对应于一个有向图的一个顶点,如果两个胞之间存在一步可达关系,则相应的顶点之间存在一条有向边.因而胞映射动力系统就完全等价于一个有向图.进而,对一个胞映射动力系统的分析就变为对一个有向图的分析.在对有向图进行分析时,可以得到图的强连通子图,强连通子图是指有向图顶点集的一个子集合,该子集合中的任何两个顶点相互可达.强连通子图对应的是相应动力系统的稳定和不稳定解集,即系统的吸引子和鞍.而不稳定解集(鞍)是其他方法不易得到的.当有向图的某一顶点不属于强连通子图时,则其对应的胞称为瞬态胞.如果一个瞬态胞可以达到某一永久自循环胞(闭的强连通子图)时,则称这个永久自循环胞为这个瞬胞的一个驻处,一个瞬态胞可以有多个驻处,按照瞬态胞所具有的驻处的数目,瞬胞被分为单驻处瞬态胞和多驻处瞬态胞,仅具有一个驻处的瞬态胞称为单驻处瞬态胞,具有多于一个驻处的瞬态胞称为多驻处瞬态胞,单驻处瞬态胞对应于吸引子的吸引域,多驻处瞬态胞对应于吸引域的边界.在Hsu[11]的论文中,详细讨论了广义胞映射,有向图和偏序集的拓扑对应,给出了实现的算法和步骤.基于Hsu[11]的工作,Xu和Hong提出了广义胞映射图方法的新版本,采用warshall算法.并将这一方法成功地用于激变的研究[12 16].本文的算法仍采用文献[11]第9章所提的算法.此外,Tongue和Gu[17]提出插值胞映射(IC M),Guder和Kreuzer[10]提出GC M的自适应方法和自适应细化方法.Sun[18,19]发表系列论文,基于短时高斯逼近研究了非线性随机系统的动力学行为,但未涉及随机分岔及不稳定解的研究.相关的研究见文献[20 25].3 参激和外激共同作用下的Duffing-van der Pol振子的随机分岔考虑如下参激和外激共同作用下的Duffing-van der Pol系统:x- x+x+x3+x2 x= 1 1(t)x+ 2 2(t),(1)式中 1, 2是两个正态白噪声过程, 1, 2为噪声的强度系数.取感兴趣的区域为D={-2 x 2,-2x 2},在这一区域内的状态胞为正规胞,区域以外的所有胞看作一个胞,称为陷胞.把区域分为200 200+1个正规胞.对于每个胞,取15 15个均匀分布的内点.对于每个内点产生50个随机样本,运用6阶龙格-库塔方法从该内点出发积分,求每个点所映射到的胞,从而在庞加莱截面上确定每个胞的象胞,则对每个胞有15 15 50=11250个轨线出发来确定其象胞.噪声的生成参见文献[26].当 1=0.0, 2= 0 0时,即不存在激励时,由图1可见,除去陷胞的影响,系统(1)有两个永久自循环胞集(吸引子);一个瞬态自循环胞集(鞍);两个单驻处瞬胞集(吸引1366物 理 学 报52卷域)和一个多驻处瞬胞集(边界).图1 系统在 1=0.0, 2=0.0时的吸引子、吸引域和鞍 图1至图12中,吸引子 的标记符号为 ;它的吸引域标记为空白;吸引子 的标记符号为 ;它的吸引域标记为空白;鞍的标记符号为 ;吸引域的边界用符号 表示考虑参数 1, 2的变化产生的随机分岔现象,讨论系统随机分岔的生成和演化,图2给出了 1- 2平面上随机分岔区域.图2 1- 2平面上的随机分岔区域计算表明,当参数 1, 2变化,在一定的参数范围内存在两次随机分岔.在区域 ,系统存在两个随机吸引子,一个随机鞍,两个随机吸引域和一个边界.在区域 ,存在一个随机吸引子,一个随机鞍,一个随机吸引域.区域 和区域 的交界形成第一次随机分岔曲线.在区域 ,存在一个随机吸引子,一个随机吸引域.区域 和区域 的交界形成第二次分岔曲线.下面分三种情况展开详细讨论.下述的讨论,为简单起见,简称随机吸引子,随机鞍,随机吸引域为吸引子,鞍,吸引域.图3 系统在 1=0.070, 2=0.000时的吸引子、吸引域和鞍图4 系统在 1=0.097, 2=0.000时的吸引子、吸引域和鞍3 1 1 0 0, 2=0.0时随机分岔的全局分析考虑仅存在参激时, 1的变化对随机系统全局特性产生的影响.随着系统参数 1的增大,吸引子及其吸引域的边界上的鞍变大,图3给出当 1=0 070时的全局图.随着系统参数 1的逐渐增大,其中的一个吸引子逐渐伸向边界上的鞍,见图4( 1=0.097),这时吸引子的指端即将接触边界上的鞍点.当随机参数 1从0 097变到0 098时,这个吸引子与在吸引域边界上的鞍发生碰撞,导致这个吸引子连同它的吸引域突然消失,在相空间原吸引子的位13676期徐 伟等:Duffing -van der Pol 振子随机分岔的全局分析置留下一个鞍,见图5( 1=0.098),此时,在我们关心的区域,仅有一个吸引子,也就是说,当系统参数 1从0 097变到0 098时,系统发生首次随机分岔,吸引子的数目由两个变为一个,吸引子的形状发生了突然变化,其中的一个吸引子变为鞍,即稳定集变为不稳定集.图5 系统在 1=0.098, 2=0.000时的吸引子、吸引域和鞍图6 系统在 1=0.099, 2=0.000时的吸引子、吸引域随着激励参数 1的进一步增大,鞍逐渐变大,并向另一个吸引子伸展,当激励参数 1由0 098变到0 099时,这个吸引子与鞍相碰撞,合而为一,突然形成一个大的吸引子,见图6( 1=0.099),即吸引子的形状突然发生变化,我们认为发生了随机分岔,此时,鞍变为吸引子,即不稳定集变为稳定集.表1给出了随参数 1变化,吸引子,吸引域,鞍的胞的变化情况.表1 随参数 1的变化吸引子、吸引域和鞍的胞数目的变化( 2=0)参数 1吸引子吸引域 吸引子吸引域 鞍0 0002516347122299020 0103016133312291040 0207015943602275250 070507145595032163080 080637141166412163090 0908001338378820883150 097940113149217598400 098961380210 00 010180 0992012379880 00 00 0图7 系统在 1=0.000, 2=0.034时的吸引子、吸引域和鞍3 2 1=0 0, 2 0.0时随机分岔的全局分析考虑仅存在外激时, 2的变化对随机系统全局特性产生的影响.随着参数 2的增大,吸引子及其吸引域的边界上的鞍变大,图7给出当 2=0 034时的全局图.随着参数 2的逐渐增大,其中的一个吸引子逐渐伸向边界上的鞍.当参数 2从0 034变到0 035时,这个吸引子与在吸引域边界上的鞍发生碰撞,导致这个吸引子连同它的吸引域突然消失,在相空间原吸引子的位置留下一个鞍,见图8( 2=0.035),此时,在我们关心的区域,仅有一个吸引子,也就是说,当系统随机参数 2从0 034变到0 035时,系统发生首次随机分岔,吸引子的数目由两个变为一个,吸引子的形状发生了突然变化,其中的一个吸引子变为鞍,即稳定集变为不稳定集.随着随机激励参数 2的进一步增大,鞍逐渐变大,并向另一个吸引子伸展,当参数 2由0 075变到0 076时,这个吸引子与鞍相碰撞,合而为一,突然形成一个大的吸引子,见图9( 2=0.075),图101368物 理 学 报52卷图8 系统在 1=0.000, 2=0.035时的吸引子、吸引域和鞍( 2=0.076),即吸引子的形状突然发生变化,此时,鞍变为吸引子,即不稳定集变为稳定集.图9 系统在 1=0.000, 2=0.075时的吸引子、吸引域和鞍图10 系统在 1=0.000, 2=0.076时的吸引子、吸引域表2给出了随参数 2变化,吸引子,吸引域,鞍的胞的变化情况.表2 随参数 2的变化吸引子、吸引域和鞍的胞数目的变化( 1=0)参数 2吸引子吸引域 吸引子 吸引域 鞍0 0338012579108218752390 03410511566117218412530 03512239257006210 03612539237006380 03713039213006570 075434382300013360 0762384376160000 0772423375770000 078245737543图11 系统在 1=0.099, 2=0.010时的吸引子、吸引域和鞍3 3 1 0 0, 2 0.0时随机分岔的全局分析考虑参激和外激共同存在时, 1, 2的变化对随机系统全局特性产生的影响.计算表明,当参数 1, 2均不为零时,类似上述的讨论,在一定的参数范围内存在两次随机分岔.但情况更为复杂一些,参激和外激的交互作用,可能使部分稳定集变为不稳定集,随着参数的增加,又可能使不稳定集变为稳定集.例如,在分岔图图2上取定 1=0.099,由图6可知,在 1=0.099, 2=0.000时,存在一个吸引子,一个吸引域.考虑参数 2的变化带来的影响,随着参数 2的增加,从吸引子中分离出一部分变为鞍,13696期徐 伟等:Duffing -van der Pol 振子随机分岔的全局分析部分稳定集变为不稳定集,这时一个吸引子,一个吸引域变为一个吸引子,一个吸引域和一个鞍.噪声起着使系统不稳定的作用,图11给出在 1=0.099, 2=0.010时的情况.随着参数 2的进一步增加,当参数 2由0 056变到0 057时,这个吸引子与鞍相碰撞,合而为一,突然形成一个大的吸引子,见图12( 1=0.099, 2=0.056),图13( 1=0.099, 2=0 057),这时部分不稳定集再次变为稳定集,又从一个吸引子,一个吸引域和一个鞍变为一个吸引子,一个吸引域,噪声起着使系统稳定的作用.表3给出了随参数 2变化,吸引子,吸引域,鞍的胞的变化情况.表3 随参数 2的变化吸引子、吸引域和鞍的胞数目的变化( 1=0 099)( 1, 2)吸引子 吸引域 吸引子吸引域鞍(0 099,0 000)201237988000(0 099,0 010)80438081001115(0 099,0 046)84437520001636(0 099,0 047)84937502001649(0 099,0 053)87737380001743(0 099,0 055)89937324001777(0 099,0 056)91037293001797(0 099,0 057)291337087000(0 099,0 058)294337057图12 系统在 1=0.099, 2=0.056时的吸引子、吸引域和鞍图13 系统在 1=0.099, 2=0.057时的吸引子、吸引域4 结 论本文引入了新的随机分岔的定义,以系统的随机吸引子(包括随机鞍)的形态(包括大小、尺寸、周期等)的突然变化表征随机系统的拓扑性质变化,用来描述随机系统的分岔.证实了对于参激和外激共同作用下的Duffing -van der Pol 系统,随机分岔主要是由于吸引子与鞍相碰撞产生的.当参数通过临界值时,这个吸引子连同它的吸引域突然消失,在原来这个吸引子的位置上留下了一个鞍,这个稳定集变为不稳定集.当随机激励参数进一步增大,通过另一个临界值时,另一个吸引子与这个鞍碰撞,鞍突然消失,这个吸引子突然变大,不稳定集变为稳定集.系统存在着两次随机分岔,图2给出了 1- 2平面上随机分岔的全局分析图.由于参激和外激交互作用,随机分岔也可能是由于当参数通过临界值时,吸引子突然分离出鞍产生的.例如,在 1=0.099, 2=0 000时,系统仅存在一个吸引子,一个吸引域,如果参数 2增加一个微小量,即 2变得不为零,则会突然从吸引子中分离出一部分变为鞍.1370物 理 学 报52卷[1]Arnold L 1998Random Dynamical Syste ms (Berli n,Berlin Heide-l berg Ne w York:Springer)1[2]Baxendale P 1986In K.Ito and T.Hida,editors Stoc hastic proce sses and their applications 1203o f LN in M athe matics ,1.[3]Crauel H and Flandoli F 1998Journal o f Dynamics and Di fferential Equations 10259[4]M eunuer C and Verga A D 1988Journal o f Statist ic al Physics 50345[5]Wol f A,Swift J B,Swinney H L,Vastano J A 1985Physica D 16285[6]Holmes P and Rand D 1980International Journal o f Non -linear Me -chanic s 15449[7]Guckenhei mer J and Holmes P 1983Nonline ar Oscillation Dynamical Syste ms and Bi furcation o f Vec to r Fields (Berlin:Springer -Verlag)[8]Namachchivaya N S 1990Applied Mathe matics and Co mputation 38101[9]Schenk -Hoppe K R 1996Nonlinear D ynamics 11255[10]Guder R and Kreuzer E 1997International Journal o f Bi furcation and Chaos .72487[11]Hsu C S 1995Inte rnational Journal o f Bi furcation and Chaos .51085[12]Xu J X and Hong L 1999Acta Me chanica Sinica 31724(in Chinese)[徐键学、洪 灵1999力学学报31724][13]Hong L and Xu J X 1999Phys .Lett .A 262361[14]Hong L and X u J X 2000Ac ta Phys .Sin .491228(i n Chi nese)[洪 灵、徐键学2000物理学报491228][15]Hong L and Xu J X 2001Acta Phys .Sin .50612(in Chinese)[洪 灵、徐键学2001物理学报50612][16]Hong L and X u J X 2002Ac ta Mec hanica Sinica 34136(in Chinese)[洪 灵、徐键学2002力学学报491228][17]Tongue B H and Gu K O 1988Journal o f Sound and Vibration 125169[18]Sun J Q and Hsu C S 1990Journal o f applied me chanics 571018[19]Sun J Q 1995Journal o f Sound and Vibration 180785[20]To C W and Li D M 1999Journal o f Sound and Vibration 219359[21]Leung H K 1998Physic al A 254146[22]He D H and Xu J X 2000Acta Phys .Sin .49833(in Chinese)[何岱海、徐键学2000物理学报49833][23]Huang Xiangao and Xu J X 2001Chin .Phys .101113[24]Li Z and Han C S 2002Chin .Phys .119[25]Lou X Z 2001Chin .Phys .1017[26]Zhu W Q,Lu M Q and Wu Q T 1993Journal o f Sound and V ib ration 165285G lobal analysis of stochastic bifu rcation in a Duffing -van der Pol system *Xu Wei 1) He Qun 1)2) Rong Ha -i Wu 3) Fang Tong 4)1)(De pa rtmen t o f Applie d Ma the ma tics ,North we ste rn Polyte chn ica l U ni versi ty ,Xi an 710072,Ch ina )2)(Department o f Appl ie d Math ematic s ,Eng inee rin g Colle ge o f Armed Polic e Force ,Xi a n 710086,China )3)(De pa rtmen t o f Applied Ma the ma tics ,Foshan Institute o f Sc ienc e an d Tech nolog y ,Foshan 528000,Ch ina )4)(In stitu te o f V ib rat ion al En gine ering ,North western Polylechn ica l Un ive rsity ,Xi an 710072,China )(Received 20September 2002;revised manu scrip t received 18November 2002)AbstractStochastic bifurcation of the Duffing -van de r Pol oscillators under both additive and multiplicative random e xcitations is stud -ied in de tail by the generalized cell mapping method using digraph.As an alternative definition,stochastic bifurcation may be de -fined as a sudde n change in charac ter of a stochastic attrac tor when the bifurcation pa ra meter of the system passe s through a c rit-i cal value.It is found that under cer tain condit ions stoc hastic bifurca tion mostly occurs when a stochastic attractor c ollides with a stochastic saddle.Our study reveals that the generalized cell mapping method with digraph is also a po we rful tool for global analy -sis of stochastic bifurcation.By this global analysis ,the mechanism of development,occurrence,and e volut ion of a stoc hastic b-i furcation can be explored clearly and vividly.Keywords :stoc hastic bifurc ation,global analysis,generalized cell mapping,stochastic attractor,stochastic saddle PACC :0547,0545*Project supported by the Nati onal Natural Science Foundati on of China (Grant No.10072049).13716期徐 伟等:Duffing -van der Pol 振子随机分岔的全局分析。
绘制Duffing振子的分叉图的程序
绘制Duffing振子的分叉图的程序这些程序思想有些可能不正确，有问题，自己改进，我不再负责对这些程序解释。
因为我都不知道道理在哪里。
但是期望您能在程序的提示下，进一步的做改进或者改正，以期获得更为精确的结果。
别照搬和迷恋别人的程序！% % %%%% 绘制Duffing振子的庞加莱截面图的程序% % buchang:已知激励下步长数值的大小, % % tend程序仿真达到150个激励周期的总时间,% clear;clc% global m c k1 k3 F0 omega%%m=1;c=0.1;k1=0;k3=1;omega=1;F0=12 % x0=[3;4];%tstart=0;Tbushu=600;buchang=(2*pi/o mega)/Tbushu;tend=(2*pi/omega)*150; % tspan=[tstart:buchang:tend];% [t,y]=ode45('dafin3',tspan,x0);%count=find(t>(2*pi/omega*40));% 去掉前40个周期的激励时间以消除瞬态响应的影响% Y=y(count,:);%TData=Y(1:Tbushu,1)-Y((Tbushu+1):Tbu shu*2,1);% [maxvalue,indices]=max(abs(TData)) %pointnumber=round((tend-2*pi/omega* 40)/buchang/Tbushu)-1;% dis=zeros(pointnumber,1);% velo=zeros(pointnumber,1);% for i=1:pointnumber%dis(i,1)=Y(Tbushu*(i-1)+indices,1);%velo(i,1)=Y(Tbushu*(i-1)+indices,2);% end% figure,plot(dis,velo,'b.','markersize',5);% %%%% 绘制Duffing振子的分叉图的程序% clear;clc% global m c k1 k3 F0 omega;% m=1;k1=0;k3=1;omega=1;F0=12;% range=[0.01:0.01:1];% YY=[];k=0;% for c=range% k=k+1;% y0=[3,4];% tspan=[0:0.01:200];% [t,Y]=ode45('dafin3',tspan,y0);% count=find(t>100);% Y=Y(count,:);% % 画x的分岔图。
duffing振子检测信号幅值的matlab代码 -回复
duffing振子检测信号幅值的matlab代码-回复Matlab代码: 振子检测信号幅值的Duffing振子模型Duffing振子是一种非线性振动系统，其运动方程可以用以下形式表示：x'' + δx' + αx + βx^3 = γcos(ωt)其中，x(t)表示振子的位移，t表示时间，α、β、γ、δ和ω是常数。
为了计算振子的振幅，我们需要首先解决振子的运动方程。
在Matlab中，可以使用ode45函数来求解常微分方程。
然后，我们可以通过对振子位移的响应进行分析，计算振幅。
以下是一步一步回答的文章，介绍如何使用Matlab代码计算Duffing振子的振幅。
第一步：定义Duffing振子的参数在编写Matlab代码之前，首先需要定义Duffing振子的参数。
这些参数包括α、β、γ、δ和ω。
可以根据具体情况选择适当的参数值。
例如，假设我们选择以下参数值：α= -1β= 1γ= 0.35δ= 0.1ω= 1.4以上参数值可用于编写Matlab代码。
第二步：编写Duffing振子的运动方程在Matlab中，可以使用函数句柄来定义Duffing振子的运动方程。
首先，我们需要定义一个函数，该函数接受振子位移和时间作为输入，并返回振子位移的导数。
以下是定义Duffing振子运动方程的Matlab代码：function dx = duffing_equation(t,x)α= -1;β= 1;γ= 0.35;δ= 0.1;ω= 1.4;dx = [x(2); -δ*x(2) - α*x(1) - β*x(1)^3 + γ*cos(ω*t)];end在上述代码中，dx表示振子位移的导数，即速度。
变量x是振子的状态变量，x(1)表示位移，x(2)表示速度。
函数返回一个包含位移和速度的列向量。
第三步：求解Duffing振子的运动方程编写完Duffing振子的运动方程后，我们可以使用ode45函数来求解振子的运动。
基于小波分解和Duffing振子的变尺度微弱信号检测
基于小波分解和Duffing振子的变尺度微弱信号检测程凯;董雪【摘要】传统的时频分析方法在对周期性微弱信号进行检测时,提取的信息具有信噪比不高的缺点,从而影响了检测效果,为此,利用Duffing振子混沌系统对噪声的强免疫力的特征,提出了一种基于小波分解和混沌阵子的混合微弱信号检测方法;首先,采用小波变换对信号进行分解,通过小波变换的平滑作用实现对含噪微弱信号的离散处理,并设计了一种根据阈值来确定分解层数的方法,然后将降噪后的重构信号作为Duffing阵子的周期驱动力并入混沌系统,采用混沌Duffing阵子阵列实现在强噪声背景下的微弱信号检测,并提出了一种临界状态策动力幅值和初始相位的自适应确定方法;在Matlab7仿真环境下进行实验,结果表明:文中方法能有效地对湮没在强噪声下的微弱信息进行检测,具有信号检测信噪比高,重构信号频率较其它方法更接近于真实频率,具有较强的可行性.【期刊名称】《计算机测量与控制》【年(卷),期】2014(022)006【总页数】3页(P1732-1734)【关键词】小波分解;微弱信号检测;噪声;混沌系统【作者】程凯;董雪【作者单位】河南教育学院信息技术系,郑州 450046;河南教育学院院长办公室,郑州 450046【正文语种】中文【中图分类】TP3120 引言从具有强噪声背景的信号中提取出有用信息，提高信噪比和抑制噪声已经成为了信号检测的重点问题［1－3］。
传统的信号处理方法往往采用重复测量取平均值的方法［4］，但其在去除噪声的同时丢失了信号的瞬态分量，因此，不能准确地分辨出信号的高频细节信息和噪声。
时频分析是处理非平稳信号的重要方法，可以提取出观测信号的时间－频率联合特征以及时域频域信息和信号频率随着时间轴变化的规律［5－6］，其缺点是只有当基函数与信号尺度函数相匹配时，才能检测出信号轮廓。
文献［7］提出了一种结合小波［8］核函数和SVM 的信号检测方法，采用谐波小波函数对信号进行降噪分析，并通过SVM 回归从而实现小样本情况下的微弱信号检测。
分岔图做法[1]
![分岔图做法[1]](https://img.taocdn.com/s3/m/ff8b4155be23482fb4da4c67.png)
>>混沌研究总结篇------一、分岔图（1.Chen系统）先打个提纲，这几天把自己混沌相关知识研究学习内容总结一下。
首先简绍几个基本概念：一、自治系统一个n阶自治的连续动态系统可以表示为可以理解为对于自治的连续系统，上相量场f是不依赖于时间t的。
二、非自治系统一个n阶非自治的连续动态系统可以表示为可以理解为对于非自治的连续系统，向量场f不仅依赖于状态变量x，而且依赖于时间t，如Duffing振子。
三、庞加莱映射庞加莱映射是一个传统的用来离散化连续系统的方法。
庞加莱映射可以用(n-1)阶的离散映射来取代n阶的连续系统。
庞加莱映射的用处正在于减小系统的阶数，并且在连续系统和离散系统之间建立了一座桥梁。
对于n阶自治系统，其对应的解对就着轨迹。
当选择作为一个(n-1)维的超平面，这样轨迹将穿越超平面。
难点主要是超平面的选取，使其对应的解穿越超平面，就可以得到一个领域内的庞加莱映射。
对于n阶非自治系统，若其外加强迫力的最小周期是T,j最终的庞加莱映射可以定义为相应的轨道P(xk)是对某个轨迹每隔T时刻采样一次获得，这种操作和每隔T时刻的频闪观测仪的行为很相似。
所以要想得到一个系统的庞加莱映射，这段话一定要好好理解，当真真知道这中间说的含义，庞加莱映射这么画其实也已经知道国。
四、分岔图分岔图的横坐标是一个变化的参数，纵坐标是你要求的某一个量的随着各参数的变化情况，而poincare则是我们选取横坐标上的某参数的某一个具体值时截面图，只不过poincare截面的选取其实可以是任意的。
下面主要研究的混沌系统有：Logistic、Henon、Lorenz、Duffing、Rossler、Chen、混沌电机模型等系统1.Chen系统先说Chen系统，因为和课题有一定的关系，而且自己以后起家也得从Chen 系统入手。
系统方程如下：dx/dt=a*(y-x)dy/dt=(c-a)*x+c*y-x*zdz/dt=x*y-b*z就是对此方程中不同参数a、b、c下对系统画分岔图，研究混沌系统(1)给定a、c，画b关于系统的分岔图结果如下图所示CODE：function fenchatuchenclc;clearXA=35;XC=28;Z=[];for XB=linspace(2,5.5,100);options = odeset('RelTol',1e-6,'AbsTol',[1e-4 1e-4 1e-5]);[T,X]=ode45('chen',[0,50],[-5 0 5],options,XA,XB,XC);n=length(X);for k=round(n/2):nif abs(X(k,1))<1Z=[Z,XB+abs(X(k,2))*i];endendendfigureplot(Z,'.','markersize',1)title('chen映射分岔图')xlabel('b'),ylabel('|x| where x=0')这组代码不完全是自己的,现在见解其中一些方法在进行自己系统的绘制，这个程序的具体原理我会在后面给出来的。
电牵引采煤机在线监测与故障诊断系统研究
2581 前言随着绿色能源风能、水能、太阳能等能源行业的进步，煤炭资源的消耗量日益减小，但由于绿色能源处于初级阶段使得煤炭资源的消耗仍十分巨大。
电牵采煤机作为煤矿生产的重要设备，其工作性能直接关系到我国煤矿的生产效率，由于矿山地质环境较为复杂，在工作过程中极易发生故障，故障如不能及时进行处理会大大降低煤矿的效率，因此对电牵引采煤机的故障诊断成为了一个重要的课题[1，2]。
此前杨成喜[3]分析了电牵引采煤机液压调高系统常见故障及其相应的解决方案。
提出提升调高系统性能的方案，为矿山提升了经济效益。
姜海[4]设计了电牵引采煤机监测系统，设计后采煤机数据记录仪实时监测采煤机的工况。
实现采煤机状态实时监测，为采煤机故障的实时监测作出一定的贡献。
本文以MG900/2215-GWD电牵引采煤机为研究对象，基于混沌和小波技术对采煤机进行振动故障分析及诊断，为矿山实验智能化开采提供一定的贡献。
2 混沌Duffing 振子在故障诊断中的应用随着机械化的不断进步，矿山采煤机的数量不断增加，同时采煤机的结构也越来越复杂，对采煤机故障进行及时的诊断也就变得十分重要。
一般来说常见的采煤机故障诊断技术主要为油液分析法、温度检测法、电气量分析法、状态参数法等。
各种分析方法均有其一定的优劣势，所以在进行采煤机故障诊断时需要考虑实际情况，合理选择适当的故障分析技术。
考虑到本文对采煤机的截割部位进行故障分析，所以选定利用振动信号对采煤机故障进行实时的监测和诊断。
振动信号监测不仅可以对部件故障进行及时的分类，同时可以给出发生故障的类型及其故障源，从而实现实时报警。
振动信号分析法根据其分析对象可分为幅值域分析法、频域分析法和解调分析法等，随着振动分析法的发展，振动信号的诊断技术已经发展到时频分析法，考虑到对振动故障即运行过程中的不平稳阶段进行分析，所以选定时频分析法进行分析，对Duffing混沌振子监测振动信号进行研究。
在实际煤矿运行过程中采煤机的振动频率是很难测量的，所以需要利用Duffing振子对振幅进行监测，振子的检测振幅及频率流程如图1所示。
二进制调制信号Duffing振子检测方法及其抗噪声机理分析
摘
要 ： 据二 进制 调制 信 号 的特 点及 Ｄ ｆｎ 振 子微 弱信 号检 测技 术 ， 究 了二 进 制 调 制 信 号 的 根 ｕ ｇ ｉ ｆ 研
周期与色噪声联合作用下分数阶Duffing_振子非平稳响应的无记忆方法
第 36 卷第 4 期2023 年8 月振 动 工 程 学 报Journal of Vibration EngineeringVol. 36 No. 4Aug. 2023周期与色噪声联合作用下分数阶Duffing振子非平稳响应的无记忆方法李书进1，张志聪1，孔凡2，韩仁杰1（1.武汉理工大学土木工程与建筑学院，湖北武汉 430070；2.合肥工业大学土木与水利工程学院，安徽合肥 230009）摘要: 基于统计线性化提出了一种求解周期与色噪声激励联合作用下分数阶Duffing系统非平稳响应的无记忆方法。
将系统响应分解为确定性周期和零均值随机分量之和，则原非线性运动方程可等效地化为一组耦合的、分别以确定性和随机动力响应为未知量的分数阶微分方程。
利用无记忆化方法将确定性和随机分数阶微分方程转化为相应的常微分方程。
利用统计线性化方法处理随机常微分方程，得到关于随机响应二阶矩的李雅普诺夫方程。
利用数值算法联立求解李雅普诺夫微分方程和确定性常微分方程。
通过Monte Carlo模拟，验证此方法的适用性和精度。
关键词: 非平稳响应；分数阶系统；无记忆方法；统计线性化；联合激励中图分类号: O324； O322 文献标志码: A 文章编号: 1004-4523（2023）04-0923-11DOI：10.16385/ki.issn.1004-4523.2023.04.005引言分数阶微积分近几十年在工程界得到广泛的关注［1］。
分数阶微积分的一个重要工程应用是黏弹性材料的力学模型的建立。
与标准线性固体模型（Standard Linear Solid Model）相比，分数阶导数模型能以较少的参数拟合实验获得黏弹性松驰和蠕变数据［2］；此外，分数阶微分方程可很好地描述动力激励下装备黏弹性控制装置（如：天然橡胶支座［3］、黏滞阻尼器［4］和黏弹性阻尼器［5⁃6］等）的动力行为。
因此，研究这类分数阶运动方程的求解方法成为学者们关注的重点。
duffing振子检测信号幅值的matlab代码
duffing振子检测信号幅值的matlab代码一、引言Duffing振子是一种广泛应用于物理、力学和工程领域的振动系统，其动态特性在许多实际应用中具有重要意义。
为了更好地理解和控制Duffing振子的行为，对其检测信号幅值进行分析是非常必要的。
本文将介绍一种基于Matlab的Duffing振子检测信号幅值的方法，并给出相应的代码实现。
二、算法介绍1. 数据采集：首先，需要使用Matlab的传感器或控制器接口采集Duffing振子的振动数据，通常包括时间序列和对应的振动幅度。
2. 信号处理：对采集到的振动数据进行预处理，包括去除噪声、平滑化等操作，以便更好地分析信号幅值。
3. 幅值计算：通过设定阈值或其他算法，计算每个时间点的振动幅度，并将其存储在新的数据结构中。
三、Matlab代码实现以下是一个简单的Matlab代码示例，用于检测Duffing振子的信号幅值：```matlab% 读取振动数据文件vibration_data = readtable('vibration_data.csv');% 设定阈值，用于计算信号幅值threshold = 0.5;% 循环遍历每个时间点for i = 1:vibration_data.size% 获取当前时间点和对应的振动幅度time = vibration_data(i, 1);amplitude = vibration_data(i, 2);% 判断振动幅度是否超过阈值，并记录结果if amplitude > thresholddetected_amplitude = amplitude;fprintf('Detected amplitude at time %f: %f\n', time, detected_amplitude);elsefprintf('No detected amplitude at time %f\n', time);endend```四、代码解析与优化这段代码实现了基本的Duffing振子检测信号幅值的功能，但在实际应用中可能需要进行一些优化和改进：1. 数据格式：确保振动数据文件格式正确，可以使用Matlab内置的读取函数或第三方工具箱来读取数据。
非线性单摆的研究
摘要单摆是日常生活中常见的一种物理现象，用一根细绳的一端拴着一个重物，把另一端固定，当重物来回摆动时，就形成了一个单摆模型。
本文讨论理想单摆和非线性单摆的分析方法，着重讨论非线性单摆的角度和角速度的关系及用摄动法求解一类特殊非线性单摆（duffing振子）。
并介绍几种常用的数学求解单摆方程的方法。
关键字：无阻尼；周期强迫；任意角度；阻尼振子；非线性；摄动法；平均法-IAbstractSingle pendulum is a common physic phenomenon in our life. Tie an object with a line, fasten the other side. When the object rock around, we can get a single pendulum.The main content of this paper is to discuss the method of analysis the ideal single pendulum and the nonlinear single pendulum. We put our eye on the relationship between the angle speed and the angle acceleration. Then we will use perturbation method to solve one special angle pendulum equation. At last we will introduce some common ways to solve the problem.Keywords:no damp, period forced, any angle, damp flap, nonlinear,perturbation method, average method.-II目录一、无阻尼振荡的分析 (1)二、周期强迫振动的分析 (4)三、摆角为任意角度的分析 (5)四、阻尼振子的分析 (8)五、有摩擦强迫振动的分析 (10)六、非线性振子的分析 (12)七、摄动法求解duffing振子方程（perturbation method） (15)1、正规摄动法（regular perturbation method） (16)2、Poicarè法： (17)八、用平均法求解单摆方程 (19)参考文献 (21)附件 (22)-III-1一、无阻尼振荡的分析如图所示，忽略细绳重量，也不计小球受到的空气阻力，则上诉单摆可看成理想单摆，对其进行受力分由牛顿第二定律得：θsin mg ma -=（1）因为2222dtd l dt s d a θ== )(θl s = （2） 把（2）代入（1）式可得 0sin 22=+θθmg dtd ml （3）将（3）两端同除以ml 可得 0sin 22=+θθlgdtd （4）令lg=0ω，其中0ω为自然频率.则（4）可变为 0sin 2022=+θωθdtd （5） 当θ很小时,θθ=sin 故有,02022=+θωθdt d （6）解此方程得:ti ti eC e C t 0021)(ωωθ-+= （7）若θ为实数,则有θθ=*,即t i t i ti t i e C e C e C e C 000021*2*1ωωωω--+=+ （8）所以, *21C C =, *12C C = （9）令ϕi e A C 21=,ϕi e AC -=22. 则有())cos(2)(0)()(0ϕωθϕωϕω+=+=+-+t A e e A t t i t i （10）-2)cos()(0ϕωθ+=t A t （11）从能量守恒方面考虑:0022=+θωθdt d 可变形为 020=+⋅⎪⎭⎫⎝⎛θωθθθdtd d dt d d （12） 令dtd θθ='，则有 0''20=+θωθθθd d （13）两边同时乘以θd ，得到 0''20=+θθωθθd d （14）在对两边求积分， ⎰⎰⎰=+θθθωθθd d d 0''20 （15）积分结果为 E =+220221'21θωθ （16）令2'21θ=T (动能),22021θω=V (势能). 则有E V T =+,机械能守恒.E =+220221'21θωθ为椭圆方程：00.20.40.60.811.2 1.4 1.6 1.82tt h e t a /d t h e t a图1 a-3xy图 1 b 图1 a 的图像为摆动角度θ及角度θ的导数随时间变化的曲线，其中实线表示角度θ随时间的变化，点线表示θ的导数（即角速度）随时间的变化。
基于Duffing振子的微弱信号检测方法及其电路的实现_赵文礼
2
Duffing 振子动力学行为
非自治的 Duffing 振子方程可表示为： d 2x dx k x x 3 F cos t （ 1） dt dt 2 式中， k dx dt 为阻尼项， x x 3 为非线性恢复力， F cos t 为驱动力。当式（ 1 ）中阻尼项和驱动力 都为零时，方程变为：
图4 频率为 60HZ 正弦信号对系统的控制 图3 频率为 20Hz 的正弦信号对系统的控制 图2
（ 8）
Duffing 振子检测电路
第6期
赵文礼等：基于 Duffing 振子的微弱信号检测方法及其电路的实现
123
弦控制信号，关键的是找出与其混沌的阈值相对应的电压幅值。比如频率为 20Hz 时，其混沌阈值所 对应的电压幅值是 6.6V， 当信号幅值大于此值时系统就从混沌状态转入大尺度周期态； 而频率为 60Hz 时的混沌临界态控制电压则是 6.3V。
1
引言
利用混沌动力学原理从强噪声背景中提取微弱的周期性有用信号是当今的热门研究领域之一。
1994 年，冯奇研究了噪声对 Duffing 振子的影响，发现噪声仅仅在有限的时间内会使系统产生复杂运 动，但是最终会趋于规则性运动 [2]； 1995 年 Haykin 利用人工神经网络方法实现了混沌背景噪声中的 小信号提取 [3]。1996 年 Leung 利用 MPSV 方法进行了混沌通信系统中如何提取有用信号的研究 [4]。尔 后 Wang Guan-Yu 等人利用混沌测量系统实现了白噪声背景下信噪比低达 -66dB 的正弦信号的测量， 并 且利用吸引子所在流形的局部切空间的补空间关系，成功提取了谐波信号 [5,6]； 2004 年李月、杨宝俊 提出了在色噪声背景下 nV 级正弦信号、方波信号、周期脉冲信号的混沌测量方法 [7]。文献 [10,11]结合 Duffing 振子作了电路系统的仿真研究。但是，目前更多的研究还处于理论计算和软件仿真方面。本文 则在前人研究的基础上，着重开展了如何通过硬件电路实现微弱信号测量的研究，设计制作了基于 Duffing 振子的微弱信号检测电路，并将电路应用于噪声背景下微弱正弦信号检测的实验。
电力系统低频振荡监测的Duffing振子可停振动系统法
Duffing方程的分岔与混沌
Duffing方程的分岔与混沌
姜春霞;邬开俊
【期刊名称】《洛阳理工学院学报（自然科学版）》
【年(卷),期】2016(026)001
【摘要】Duffing方程是非线性振动系统中的一种具有代表性的微分方程式, 许多工程实际中的非线性振动问题都可以利用该数学方程来研究. 本文应用动力系统的分岔理论和混沌理论, 研究在五次非线性恢复力、一个激励和一个外力作用下的Duffing方程. 用理论解析法以及数值仿真的方式求出方程自由振动、强迫振动、混沌振动. 用数值仿真 (分岔图、相图、庞加莱映射图) 研究了该系统的复杂动态. 这些数值仿真展示了周期倍化序列到混沌、混沌和周期窗口的交替出现以及混沌的突然消失等动态行为.
【总页数】6页(P83-88)
【作者】姜春霞;邬开俊
【作者单位】辽宁省交通高等专科学校轨道交通工程系,辽宁沈阳110112;兰州交通大学电子与信息工程学院,甘肃兰州730070
【正文语种】中文
【中图分类】O29;TH113.1
【相关文献】
1.两个自由度Duffing系统的分岔及混沌分析 [J], 曹振邦;乐源
2.两个自由度Duffing系统的分岔及混沌分析 [J], 曹振邦;乐源;
3.Van der Pol-Duffing耦合系统的分岔与混沌控制 [J], 褚衍东;李险峰;张建刚
4.双边碰撞Duffing振子的对称性、尖点分岔与混沌 [J], 李冠强;谢建华
5.基于增量谐波平衡法的Mathieu-Duffing振子分岔及通往混沌道路分析 [J], 陈树辉;沈建和
因版权原因，仅展示原文概要，查看原文内容请购买。
怎么用matlab画分叉图,混沌------分岔图绘制不同方法的总结、比较（转）
怎么⽤matlab画分叉图,混沌------分岔图绘制不同⽅法的总结、⽐较（转）经过近期的研究发现，⽬前对于系统单参数分岔图的计算共有以下的⼏种⽅法：1)最⼤值法即对系统微分⽅程(组)进⾏求解，对求解的结果⽤getmax函数进⾏取点，并绘图。
2)Poincare截⾯法对系统参数的每⼀次取值，绘制其Poincare截⾯，进⽽得到其分岔图。
这种⽅法需要注意的是，⾃治系统的Poincare截⾯是选取⼀超平⾯，平⾯上点的分布即构成⼀Poincare截⾯，⾮⾃治系统的Poincare截⾯则是根据系统激励的频率进⾏取点并绘图。
本帖将以Lorenz系统为例，对这两种⽅法进⾏⽐较⾸先对第⼆种⽅法进⾏阐述。
编程如下(matlab)Lorenz系统：function dy = Lorenz(t,y)% Lorenz系统% 系统微分⽅程：% dx/dt = -a(x-y)% dy/dt = x(r-z)-y% dz/dt = xy-bz% a=y(4)% r=y(5)% b=y(6)dy=zeros(6,1);dy(1)=-y(4)*(y(1)-y(2));dy(2)=y(1)*(y(5)-y(3))-y(2);dy(3)=y(1)*y(2)-y(6)*y(3);dy(4)=0;dy(5)=0;dy(6)=0;随r的分岔图求解程序：——按照x=y平⾯取截⾯function Lorenz_bifur_rZ=[];for r=linspace(1,500,1000);% 舍弃前⾯迭带的结果，⽤后⾯的结果画图[T,Y]=ode45('Lorenz',[0,1],[1;1;1;16;r;4]); [T,Y]=ode45('Lorenz',[0,50],Y(length(Y),:));Y(:,1)=Y(:,2)-Y(:,1);% 对计算结果进⾏判断，如果点满⾜x=y，则取点for k=2:length(Y)f=k-1;if Y(k,1)<0if Y(f,1)>0y=Y(k,2)-Y(k,1)*(Y(f,2)-Y(k,2))/(Y(f,1)-Y(k,1)); Z=[Z r+abs(y)*i];endelseif Y(f,1)<0y=Y(k,2)-Y(k,1)*(Y(f,2)-Y(k,2))/(Y(f,1)-Y(k,1)); Z=[Z r+abs(y)*i];endendendendplot(Z,'.','markersize',1)title('Lorenz映射分岔图')xlabel('r'),ylabel('|y| where x=y')结果如图1所⽰。
Duffing方程介绍与仿真应用
非线性电路理论及应用报告• Duffing方程介绍与仿真应用姓名：马博学号：25班级：硕3022班完成时间：在非线性振动理论研究中，Duffing方程是一种具有代表性的微分方程式。
本文首先对Duffing方程进行了简单介绍，包括其类型以及根据电路的推导等；其次，本文对硬特性的Duffing方程进行了不同参数下的Mat lab仿真；最后，本文介绍了Duffing方程的微弱信号频率检测，以Holmes型Duffing方程为例进行了分析说明。
关镀词：Duffing方程；非线性；Mat la b仿真；混沌；弱信号检测1Duffing方程简介非线性振动问题的研究通常包括定性研究与定量研究。
定性研究的主要内容包括方程解的存在性、唯一性、周期性和稳定性的研究等。
著名的Duffing方程在非线性动力学系统的研究中占有重要地位。
其特点之一是在Duffing方程等号右边加上了外加强迫项，进而形成了非自治非线性系统。
正是由于系统的本征频率与外加周期强迫项的频率的相互作用，才使得该方程中蕴含着极其丰富的内容：倍周期分叉、混沌、清晰大周期等现象⑴。
(Duffing方程的准形式为：d2x ^dx /、、+ / +g(x) = /(x,t)dt~ dt其中5>0为阻尼系数，g(x)是含有三次方项的非线性函数，/(x,t)为一周期函数。
Duffing方程通常作如下分类⑵：1.假设g(x)满足超线性条件lim型十L T—OO x则称Duffing方程是超线性的；2.假设g(x)满足次线性条件lim 型=0L T—oc 牙则称Duffing方程是次线性的；3.假设g(x)满足半线性条件0 < lim inf < lim sup ^―— < -+<olxlT8 牙Ldoc x则称Duffing方程是半线性的。
若将Duffing方程规范化，有以下四种基本类型⑶:d2x f dx八彳八”-—T + ^ —+ X(t) + X (t) = j cos(t) (1-1)d2x . dx八彳八彳八,+x(t)x(t)=/cos(t)+(1-2)dxd2x+ k — + x3(t) = f cos(t) (1-3)crx f dx、勺八 c眉 + 匚- x(t) + F (t) = f cos(t) (1一4)其中k大于零，是阻尼系数，/cos（t）是系统外力。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
这些程序思想有些可能不正确，有问题，自己改进，我不再负责对这些程序解释。
因为我都不知道道理在哪里。
但是期望您能在程序的提示下，进一步的做改进或者改正，以期获得更为精确的结果。
别照搬和迷恋别人的程序！% % %%%% 绘制Duffing振子的庞加莱截面图的程序% % buchang:已知激励下步长数值的大小,% % tend程序仿真达到150个激励周期的总时间,% clear;clc% global m c k1 k3 F0 omega%% m=1;c=0.1;k1=0;k3=1;omega=1;F0=12% x0=[3;4];% tstart=0;Tbushu=600;buchang=(2*pi/omega)/Tbushu;tend=(2*pi/omega)*150;% tspan=[tstart:buchang:tend];% [t,y]=ode45('dafin3',tspan,x0);% count=find(t>(2*pi/omega*40)); % 去掉前40个周期的激励时间以消除瞬态响应的影响% Y=y(count,:);% TData=Y(1:Tbushu,1)-Y((Tbushu+1):Tbushu*2,1);% [maxvalue,indices]=max(abs(TData))% pointnumber=round((tend-2*pi/omega*40)/buchang/Tbushu)-1;% dis=zeros(pointnumber,1);% velo=zeros(pointnumber,1);% for i=1:pointnumber% dis(i,1)=Y(Tbushu*(i-1)+indices,1);% velo(i,1)=Y(Tbushu*(i-1)+indices,2);% end% figure,plot(dis,velo,'b.','markersize',5);% %%%% 绘制Duffing振子的分叉图的程序% clear;clc% global m c k1 k3 F0 omega;% m=1;k1=0;k3=1;omega=1;F0=12;% range=[0.01:0.01:1];% YY=[];k=0;% for c=range% k=k+1;% y0=[3,4];% tspan=[0:0.01:200];% [t,Y]=ode45('dafin3',tspan,y0);% count=find(t>100);% Y=Y(count,:);% % 画x的分岔图。
% j=1;% n=length(Y(:,1));% for i=2:n-1% if Y(i-1,1)+eps<Y(i,1) & Y(i,1)>Y(i+1,1)+eps %简单的取出局部最大值。
% YY(k,j)=Y(i,1); %使最大值计数个数自动增加% j=j+1;% end% end% if j>1% plot(c,YY(k,[1:j-1]),'b.','markersize',5);% end% hold on;% index(k)=j-1;% end% xlabel('c');% ylabel('x max');% title('dafin bifurcation diagram');% % % 绘制分岔图的程序% clear,clc% global m c k1 k3 F0 omega%% m=1;c=0.1;k1=0;k3=1;omega=1;F0=12;% ccanshu=0.01:0.01:1;% for k=1:100% c=ccanshu(k)% x0=[3;4];% tspan=[0:0.01*2*pi:500];% [t,y]=ode45('dafin3',tspan,x0);% dis=zeros(50,1);% velo=zeros(50,1);% for i=1:50% dis(i,1)=y(100*(i+20),1);% velo(i,1)=y(100*(i+20),2);% end% Dismatrix(k,:)=dis';% end% figure,plot(ccanshu,Dismatrix,'b.','markersize',3); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%% %线性参数k1的变化产生的分岔图% clear;clc% global m c k1 k3 F0 omega% m=1;c=0.1;k3=1;omega=1;F0=12;% kcanshu=0.01:0.01:2;% k1=kcanshu(k)% x0=[3;4];% tspan=[0:0.01*2*pi:500];% [t,y]=ode45('dafin3',tspan,x0);% dis=zeros(50,1);% velo=zeros(50,1);% for i=1:50% dis(i,1)=y(100*(i+20),1);% velo(i,1)=y(100*(i+20),2);% end% Dismatrix(k,:)=dis';% end% plot(kcanshu,Dismatrix,'b.','markersize',5); % title('参数变化下的分岔图')% xlabel('线性刚度参数k1的变化')% ylabel('X值')% %非线性参数k3的变化产生的分岔图% clear;clc% global m c k1 k3 F0 omega% m=1;c=0.1;k1=0;omega=1;F0=12;% kcanshu=0.01:0.01:2;% for k=1:200% k3=kcanshu(k)% x0=[3;4];% tspan=[0:0.01*2*pi:500];% [t,y]=ode45('dafin3',tspan,x0);% dis=zeros(50,1);% velo=zeros(50,1);% for i=1:50% dis(i,1)=y(100*(i+20),1);% velo(i,1)=y(100*(i+20),2);% end% Dismatrix(k,:)=dis';% end% plot(kcanshu,Dismatrix,'b.','markersize',5); % title('参数变化下的分岔图')% xlabel('非线性参数k3的变化')% ylabel('X值')% %激励参数F0变化产生的分岔图% clear;clc% global m c k1 k3 F0 omega% m=1;c=0.1;k1=0;k3=1;omega=1;% F0canshu=0.1:0.1:20;% F0=F0canshu(k)% x0=[3;4];% tspan=[0:0.01*2*pi:500];% [t,y]=ode45('dafin3',tspan,x0);% dis=zeros(50,1);% velo=zeros(50,1);% for i=1:50% dis(i,1)=y(100*(i+20),1);% velo(i,1)=y(100*(i+20),2);% end% Dismatrix(k,:)=dis';% end% plot(F0canshu,Dismatrix,'b.','markersize',5);% title('参数变化下的分岔图')% xlabel('激励参数F0的变化')% ylabel('X值')% % %激励频率omega变化产生的分岔图% clear;clc% global m c k1 k3 F0 omega% m=1;c=0.1;k1=0;k3=1;F0=12;% omegacanshu=0.1:0.1:10;% for k=1:100% omega=omegacanshu(k)% x0=[3;4];% tspan=[0:0.01*2*pi/omega:500];% [t,y]=ode45('dafin3',tspan,x0);% dis=zeros(50,1);% velo=zeros(50,1);% for i=1:50% dis(i,1)=y(round(100*omega*(i+20)),1); % velo(i,1)=y(round(100*omega*(i+20)),2); % end% Dismatrix(k,:)=dis';% end% plot(omegacanshu,Dismatrix,'b.','markersize',5); % title('参数变化下的分岔图')% xlabel('激励频率omega的变化')% ylabel('X值')% clear;clc% global m c k1 k3 F0 omega% n=3,rhs_ext_fcn=@dafin_ext2,fcn_integrator=@ode45,tstart=0,stept=0.5,tend=200,% ystart=[3 4 0],ioutp=10,% m=1;c=0.1;k1=0;F0=12;k3=1;% omegacanshu=0.1:0.1:10;% for k=1:100% omega = omegacanshu(1,k),lyapunovzhishu(k,:)=lyapunovfun(n,rhs_ext_fcn,fcn_integrator,tstart,stept,tend,y start,ioutp)% end% figure,plot(omegacanshu,lyapunovzhishu),% title('Lyapunov 动力学指数');% xlabel('激励频率omega变化'); ylabel('Lyapunov 指数');% % % 绘制分岔图的程序% clear;clc% global m c k1 k3 F0 omega%% m=1;c=0.1;k1=0;k3=1;omega=1;F0=12;% ccanshu=0.01:0.01:1;% for k=1:100% c=ccanshu(k)% x0=[3;4];% tstart=0;Tbushu=100;buchang=(2*pi/omega)/Tbushu;tend=(2*pi/omega)*200;% tspan=[tstart:buchang:tend];% [t,y]=ode45('dafin3',tspan,x0);% count=find(t>(2*pi/omega*40)); % 去掉前40个周期的激励时间以消除瞬态响应的影响% Y=y(count,:);% TData=Y(1:Tbushu,1)-Y((Tbushu+1):Tbushu*2,1);% if k==1% [maxvalue,indices]=max(abs(TData));% end% pointnumber=round((tend-2*pi/omega*40)/buchang/Tbushu)-1;% dis=zeros(pointnumber,1);% velo=zeros(pointnumber,1);% for i=1:pointnumber% dis(i,1)=Y(Tbushu*(i-1)+indices,1);% velo(i,1)=Y(Tbushu*(i-1)+indices,2);% end% Dismatrix(k,:)=dis';% end% plot(ccanshu,Dismatrix,'b.','markersize',3);% % % 绘制分岔图的程序% clear,clc% global m c k1 k3 F0 omega% m=1;c=0.4;k1=-1;k3=1;F0=3;omega=2;% ccanshu=0.01:0.01:1;% for k=1:100% c=ccanshu(k)% x0=[2;0];% tstart=0;Tbushu=100;buchang=(2*pi/omega)/Tbushu;tend=(2*pi/omega)*200;% tspan=[tstart:buchang:tend];% [t,y]=ode45('dafin3',tspan,x0);% count=find(t>(2*pi/omega*40)); % 去掉前40个周期的激励时间以消除瞬态响应的影响% Y=y(count,:);% TData=Y(1:Tbushu,1)-Y((Tbushu+1):Tbushu*2,1);% if k==1% [maxvalue,indices]=max(abs(TData));% end% pointnumber=round((tend-2*pi/omega*40)/buchang/Tbushu)-1;% dis=zeros(pointnumber,1);% velo=zeros(pointnumber,1);% for i=1:pointnumber% dis(i,1)=Y(Tbushu*(i-1)+indices,1);% velo(i,1)=Y(Tbushu*(i-1)+indices,2);% end% Dismatrix(k,:)=dis';% end% figure,plot(ccanshu,Dismatrix,'b.','markersize',3);% set(gca,'fontsize',20);% title('随参数变化的分岔图','fontsize',20);% xlabel('随阻尼参数c变化','fontsize',20);% % % 绘制分岔图的程序% clear,clc% global m c k1 k3 F0 omega% m=1;c=0.4;k1=-1;k3=1;F0=3;omega=2;% k1canshu=-1:0.01:0.99;% for k=1:200% k1=k1canshu(k)% x0=[2;0];% tstart=0;Tbushu=100;buchang=(2*pi/omega)/Tbushu;tend=(2*pi/omega)*200;% tspan=[tstart:buchang:tend];% [t,y]=ode45('dafin3',tspan,x0);% count=find(t>(2*pi/omega*40)); % 去掉前40个周期的激励时间以消除瞬态响应的影响% Y=y(count,:);% TData=Y(1:Tbushu,1)-Y((Tbushu+1):Tbushu*2,1);% if k==1% [maxvalue,indices]=max(abs(TData));% end% pointnumber=round((tend-2*pi/omega*40)/buchang/Tbushu)-1;% dis=zeros(pointnumber,1);% velo=zeros(pointnumber,1);% for i=1:pointnumber% dis(i,1)=Y(Tbushu*(i-1)+indices,1);% velo(i,1)=Y(Tbushu*(i-1)+indices,2);% end% Dismatrix(k,:)=dis';% end% figure,plot(k1canshu,Dismatrix,'b.','markersize',3);% axis([-1,1,-1,4])% set(gca,'fontsize',20);% title('随参数变化的分岔图','fontsize',20);% xlabel('随线性刚度参数k1的变化','fontsize',20);% % % 绘制分岔图的程序% clear,clc% global m c k1 k3 F0 omega% m=1;c=0.4;k1=-1;k3=1;F0=3;omega=2;% k3canshu=0.01:0.01:1;% for k=1:100% k3=k3canshu(k)% x0=[2;0];% tstart=0;Tbushu=100;buchang=(2*pi/omega)/Tbushu;tend=(2*pi/omega)*200;% tspan=[tstart:buchang:tend];% [t,y]=ode45('dafin3',tspan,x0);% count=find(t>(2*pi/omega*40)); % 去掉前40个周期的激励时间以消除瞬态响应的影响% Y=y(count,:);% TData=Y(1:Tbushu,1)-Y((Tbushu+1):Tbushu*2,1);% if k==1% [maxvalue,indices]=max(abs(TData));% end% pointnumber=round((tend-2*pi/omega*40)/buchang/Tbushu)-1;% dis=zeros(pointnumber,1);% velo=zeros(pointnumber,1);% for i=1:pointnumber% dis(i,1)=Y(Tbushu*(i-1)+indices,1);% velo(i,1)=Y(Tbushu*(i-1)+indices,2);% end% Dismatrix(k,:)=dis';% end% figure,plot(k3canshu,Dismatrix,'b.','markersize',3);% set(gca,'fontsize',20);% title('随参数变化的分岔图','fontsize',20);% xlabel('随非线性刚度参数k3的变化','fontsize',20);% % % 绘制分岔图的程序% clear,clc% global m c k1 k3 F0 omega% m=1;c=0.4;k1=-1;k3=1;F0=3;omega=2;% F0canshu=0.1:0.1:10;% for k=1:100% F0=F0canshu(k)% x0=[2;0];% tstart=0;Tbushu=100;buchang=(2*pi/omega)/Tbushu;tend=(2*pi/omega)*200;% tspan=[tstart:buchang:tend];% [t,y]=ode45('dafin3',tspan,x0);% count=find(t>(2*pi/omega*40)); % 去掉前40个周期的激励时间以消除瞬态响应的影响% Y=y(count,:);% TData=Y(1:Tbushu,1)-Y((Tbushu+1):Tbushu*2,1);% if k==1% [maxvalue,indices]=max(abs(TData));% end% pointnumber=round((tend-2*pi/omega*40)/buchang/Tbushu)-1;% dis=zeros(pointnumber,1);% velo=zeros(pointnumber,1);% for i=1:pointnumber% dis(i,1)=Y(Tbushu*(i-1)+indices,1);% velo(i,1)=Y(Tbushu*(i-1)+indices,2);% end% Dismatrix(k,:)=dis';% end% figure,plot(F0canshu,Dismatrix,'b.','markersize',3);% set(gca,'fontsize',20);% title('随参数变化的分岔图','fontsize',20);% xlabel('随外界激励幅值F0的变化','fontsize',20);% %激励频率omega变化产生的分岔图clear;clcglobal m c k1 k3 F0 omegam=1;c=0.1;k1=0;k3=1;F0=12;omegacanshu=0.1:0.1:10;for k=1:100omega=omegacanshu(k)x0=[3;4];tstart=0;Tbushu=100;buchang=(2*pi/omega)/Tbushu;tend=(2*pi/omega)*200;tspan=[tstart:buchang:tend];[t,y]=ode45('dafin3',tspan,x0);count=find(t>(2*pi/omega*40)); % 去掉前40个周期的激励时间以消除瞬态响应的影响Y=y(count,:);TData=Y(1:Tbushu,1)-Y((Tbushu+1):Tbushu*2,1);if k==1[maxvalue,indices]=max(abs(TData));endpointnumber=round((tend-2*pi/omega*40)/buchang/Tbushu)-1;dis=zeros(pointnumber,1);velo=zeros(pointnumber,1);for i=1:pointnumberdis(i,1)=Y(Tbushu*(i-1)+indices,1);velo(i,1)=Y(Tbushu*(i-1)+indices,2);endDismatrix(k,:)=dis';endfigure,plot(omegacanshu,Dismatrix,'b.','markersize',3);set(gca,'fontsize',20);title('随参数变化的分岔图','fontsize',20);xlabel('随外界激励频率omega的变化','fontsize',20);% %%%%%%%%%%%%%%%%%%%%%% clear,clc% n=3,rhs_ext_fcn=@dafin_ext2,fcn_integrator=@ode45,tstart=0,stept=0.5,tend=200,% ystart=[2 0 0],ioutp=10,% global m c k1 k3 F0 omega% m=1;c=0.4;k1=-1;k3=1;F0=3;omega=2;% kcanshu=-1:0.01:0.99;% for k=1:200%k1=kcanshu(1,k),lyapunovzhishu(k,:)=lyapunovfun(n,rhs_ext_fcn,fcn_integrator,tstart,stept,tend,yst art,ioutp)% end% figure,plot(kcanshu,lyapunovzhishu),% title('Lyapunov 动力学指数');% xlabel('线性刚度参数k1,角频率为2rad/s'); ylabel('Lyapunov 指数');%%% clear,clc% n=3,rhs_ext_fcn=@dafin_ext2,fcn_integrator=@ode45,tstart=0,stept=0.5,tend=200,% ystart=[2 0 0],ioutp=10,% global m c k1 k3 F0 omega% m=1;c=0.4;k1=-1;k3=1;F0=3;omega=2;% kcanshu=[0.1:0.1:10,10:10:1000];% for k=1:200%k3=kcanshu(1,k),lyapunovzhishu(k,:)=lyapunovfun(n,rhs_ext_fcn,fcn_integrator,tstart,stept,tend,yst art,ioutp)% end% figure,plot(kcanshu,lyapunovzhishu),set(gca,'fontsize',20);% title('Lyapunov 动力学指数','fontsize',20);% xlabel('非线性刚度参数k3,角频率为2rad/s','fontsize',20); ylabel('Lyapunov 指数','fontsize',20);。