自激振荡系统matlab仿真课程设计

题目 自激振动系统的MATLAB仿真专业姓名学号主要内容、基本要求、主要参考资料等 主要内容：x2丿dt所描述的非线性有阻尼研究范德波耳(Van der pol)方程d-2-dt 2的自激振动系统，其中卩是一个小的正的参量，x0是常数。下面简称范德波耳方程为VDP 方程。在VDP方程中，增加外驱动力V cos wt项所得到的方程电代2_ x 2独+ w 2 x + V cos=0称强迫V DP方程，其中外驱动力的振幅、角频率分别是dt 20dt 0V和w。试研究强迫VDP方程的行为。基本要求：1演示VDP方程所描述的系统在非线性能源供给下，从任意初始条件出发都能产生稳定 的周期性运动。2采用庞加莱映像，演示强迫VDP方程在不同参数下所存在四种吸引子，即周期1吸引 子、周期2吸引子、不变环面吸引子和奇怪吸引子。3对于强迫V DP方程，在v和w为定值条件下，逐渐增大“值，将出现周期倍分岔和混沌 现象。主要参考资料：1 Steven E. Koonin,秦克诚译.计算物理学.北京：高等教育出版社，1993.2 马文淦等.计算物理学.合肥：中国科学技术大学出版社，1992.3 张志涌.精通MATLAB6.5.北京：北京航空航天大学出版社，2003.完成期限指导教师专业负责人2012年 2 月 23 日第1 章 概述1.1 自激震荡自激振动是一种对科学技术非常有意义而在自然界又广泛存在的非线性振 动. 我们知道,对线性阻尼振动系统 ,严格的周期运动只能由受周期性驱动力作用 的受迫振动产生;而对非线性系统,有一种自激振动系统,在非振动的能源供给下, 它能产生严格的周期运动,这是人们十分感兴趣的现象。自激系统是一个非线性 的有阻尼的振动系统 ,在振动过程中伴随有能量损耗 ,但系统存在一种机制 ,使能 量能够由非振动能源通过系统本身的反馈调节,及时适量地得到补充,从而产生一 个稳定的不衰减的周期运动,这样的振动称为自激振动1。自激振动现象是一种普遍现象。如钟摆、弦乐器以及人的心脏的周期性跳动。 活塞发动机的周期性运动等都是利用这种现象来建立不衰减的周期运动；但有些 自激振动是十分有害的，这些现象应该设法避免。1.2范德波耳(Van der pol )方程在此课程设计中，我们主要研究范德波耳(Van der pol)方程:d2 x()dx(1-1)一口 2 - x2 丿+o 2 x - 0dt 20dt0所描述的非线性有阻尼的自激振动系统，其中卩是一个小的正的参量，x0是 常数。下面简称范德波耳方程为VDP方程。在VDP方程中，增加外驱动力V cos ot项所得到的方程:d 2 x ( 一卩V dt 2)dxdt1-2)称强迫VDP方程，其中外驱动力的振幅、角频率分别是V和o。VDP 方程所描述的系统在非线性能源供给下，从任意初始条件出发都能产 生稳定的周期性运动。而对于强迫VDP方程，在V和为o定值条件下，卩逐渐增大，将出现周期倍分岔和混浊现象。第 2 章 自激振荡系统的周期性运动2.1 课题解析自激系统是一个非线性有阻尼的振动系统，在运动过程中伴随有能量损耗， 但系统存在一种机制，使能量能够由非振动的能源通过系统本身的反馈 调节， 及时适量地得到补充，从而产生一个稳定的不衰减的周期运动，这样的振动称为 自激振动。对VDP方程，可从机械振动角度理解，2 - x 2）是阻尼系数，它是变化0的，如果x |x |，则阻尼系数为正，系统将受阻尼，能量将逐渐减少，但如果0x |x I，则发生负阻尼，意味着不仅不消耗系统的能量，反而给系统提供能量。0此系统能通过自动的反馈调节，使得在一个振动过程中，补充的能量正好等于消 耗的能量，从而系统作稳定的周期振动。2.2 MATLAB 程序设计及演示取方程（1-1）中的x = 1，O = 1， 0.3，0.66,，0.85，1.08这四个初始条00件进行编程，程序详见附录-程序（1）。2.3 MATLAB 演示结果和分析根据图 3.1 显示的结果,我们做出以下结论:给出任一初始条件，通过计算机 数值求解可以证明它的相轨道都将趋向于一条闭合曲线，这一条闭合曲线，成为 极限环，极限环以外的相轨道向里盘旋，而极限环以内的相轨道则向外盘旋，都 趋向极限环，说明不论初始情况如何，系统最终都到达以极限环描述的周期性运 动2。图 3.1 四种初始条件下的轨道第 3 章 强迫 VDP 方程在不同参数下的四种吸引子下面研究强迫 VDP 方程的行为，我们同时采用时间历程图、相图、庞加莱 映像图来研究系统在不同参数条件下的动力学行为，可以看到存在不同的吸引 子，即周期 1吸引子，周期2吸引子，不变环面吸引子和奇怪吸引子。3.1 庞加莱映射为了更清楚地了解运动的形态，庞加莱对连续运动的轨迹用一个截面(叫庞 加莱截面)将其横截，那么根据轨迹在截面上穿过的情况，就可以简洁地判断运 动的形态，由此所得图像叫庞加莱映像。在截面图上，轨迹下一次 穿过截面的 点 x 可以看成前一次穿过的点 x 的一种映射。n +1nx = f (X)(n=0,l,2,.)(3-1)n + 1n这个映射就叫庞加莱映射。它把一个连续的运动化为简洁的离散映射来研 究。在庞加莱映射中的不动点反映了相空间的周期运动，如果运动是二倍周期的， 则庞加莱映射是两个不动点，四倍周期则有四个不动点等。绘制庞加莱映射是在普通的相平面上进行，它不是像画相轨道那样随时间变 化连续地画出相点，而是每隔一个外激励周期(T = 2兀/)取一个点，例如取 样的时刻可以是t=0,T,2T相应的相点记为P (X , y八P (X , y )、P (X , y )这 0 0 0 1 1 1 2 2 2 些离散相点就构成了庞加莱映射3。设 y = x， y12dX，则(1-2)式可化为 dtdt3-2)dy尸=y (x2 - y2) y - w 2 y - V cos w tdt 0 1 2 0 13.2强迫VDP方程在不同参数下的种吸引子对于（3-2）式dy=卩(X 2 - y 2)y - W 2 y - V COS W tdt 0 1 2 0 1取X 2 = 1，W 2 = 1，进行以下数值计算研究：00在卩二0.85，V = 1，= 0.44条件下，存在周期1吸引子，它的周期等于外激 励的周期，代表主谐波运动，如图 3.1所示：r IJL1LIdLIDtoo20Q300400500SOO7D0600虎垃菜貳面一周期勺吸引子图3.1周期1吸引子W = 0.441.02 ， V = 1 ，它的周期等于外激励的整数倍，代表次谐波运动，如图 3.2所示：62鹿加粟載而一伺期凍引子位移曲线在卩=0.66，V = 1，co = 0.44条件下，存在不变环面吸引子，它代表准周期（拟周期）运动，如图 3.3所示图 3.3 不变环面吸引子R = 1.08，V = 1，o= 0.44条件下存在奇怪吸引子，它代表混浊运动。 如图 3.4 所示图 3.4 奇怪吸引子保持V和o为定值，逐渐增大卩，将显示系统状态演化过程全貌的图。而前 四种情况中，看到的只是卩取4个值的片断情况，图形显示，当卩由0.9连续变 化到 1.2 时，系统运动状态逐渐由周期 1 过渡到周期 2（发生了周期倍分岔）再 过渡到混浊状态。如图 3.5 所示。3.3 MATLAB 演示程序设计在程序中，这几种过程的计算是相同的，所以用for循环来完成前面四种计 算，程序详见附录-程序（2）。 计算中在每个外激励周期内计算 1000 个相点， 为了作出庞加莱映射，每隔1000个点保留一个点数据,所以程序运行的时间较长， 对第五种情况，由于计算量大，将它另外编写一个程序，程序详见附录-程序（3）。 计算中在每个周期内计算 100 个相点，庞加莱映射是每隔 100 个点保留一个点数 据。第4章 总结计算物理学中，应用MATLAB软件求解具体问题的例子多不胜数。它以编 程效率高、使用方便、扩充能力强、语句简单内涵丰富等特点广泛应用于数据分 析、数值与符号计算、图像与数字信号处理领域。应用MATLAB软件对自激振 动系统进行仿真模拟，可以直观的显示出自激振动系统的周期性运动，以及强迫 VDP方程在不同参数下所存在四种吸引子的特点，是可以观察到对于强迫VDP 方程，在v和 为定值条件下，逐渐增大值，将出现的周期倍分岔和混沌现象。通过应用 MATLAB 软件对自激振动系统进行仿真模拟，我们对自激振动 系统有了更加深刻的了解。通过演示程序的设计编程，使我们更加熟练地掌握了 MATLAB 软件的使用方法。在编程过程中，细微处符号语句的使用让我们吃够 了苦头，往往差之毫厘失之千里。此外我们也了解了庞加莱映射等相关知识，对 于计算物理学也有了更深刻的认识。附录程序(1)：u=0.3,0.66,0.85 1.08;x0=1; w0=1; v=0; w=0.44;T=2*pi/w;for j=1:4t,y=ode23(zjzdfun,0:T/1000:50*T,4,4,u(j),x0,w0,v,w); figureplot(y(3000:end,1),y(3000:end,2);axis(-3 3 -4 4) xlabel(x);ylabel(y);end程序(2)：u=0.85, 1.02, 0.66, 1.08;x0=1; w0=1; v=1; w=0.44;T=2*pi/w;strl=庞加莱截面一周期1吸引子；str2=庞加莱截面一周期2吸引子；str3=庞加莱截面一不变环面吸引子；str4=庞加莱截面一奇怪吸引子；for j=l:4t,y=ode23(zjzdfun,0:T/l000:50*T,4,4,u(j),x0,w0,v,w)； figuresubplot(2,l,l) plot(t,y(:,l)；title(位移曲线)；xlabel(x)；ylabel(v)；subplot(2,2,3)plot(y(3000:end,l),y(3000:end,2)；axis(-3 3 -4 4)xlabel(x)；ylabel(v)；title (相图)；subplot(2,2,4) axis(-3 1 -1 1) hold on for i=7000:1000:14000 plot(y(i,1),y(i,2),r.);end title(strj); end程序(3)：u=0.8:0.001:1.2; v=1; x0=1;w0=1; w=0.44;T=2*pi/w; axis(0.9 1.2 -0.8 1) hold on for j=1:length(u)t,y=ode23(zjzdfun,0:T/100:70*T,4,4,u(j),x0,w0,v,w); plot(u(j),y(500:100:1400,2),linewidth,2);end 程序(4)(该程序为调用函数)：function ydot=vdbfun(t,y,flag,u,x0,w0,v,w) ydot=y(2);u*(x0人2-y(l)人2)*y(2)-y(l)*w02-v*cos(w*t);参考文献1 丁文静.自激振动，清华大学出版社，2009，(34)：325-328.2 胡静,彭芳麟，管靖，卢圣治理论力学中非线性问题的MATLAB数值解，大学 物理COLL EGE PHYSICS 第20卷第 10期,20013 张志涌.精通MATLAB6.5北京：北京航空航天大学出版社，2003.