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

课程 计算物理和MATLAB 课程设计 题目 自激振动系统的MATLAB 仿真 专业
姓名
学号
主要内容、基本要求、主要参考资料等 主要内容：
研究范•德•波耳(Van der pol)方程()0202202
2
=+--x dt
dx x x dt x d ωμ所描述的非线性有阻尼
的自激振动系统，其中μ是一个小的正的参量，0x 是常数。

下面简称范•德•波耳方程为VDP 方程。

在VDP 方程中，增加外驱动力t V ωcos 项所得到的方程
(
)0
cos 2
02
2022
=++--t V x dt
dx
x
x dt
x d ωωμ称强迫VDP 方程，其中外驱动力的振幅、角频率分别是
V 和ω。

试研究强迫VDP 方程的行为。

基本要求：
1.演示VDP 方程所描述的系统在非线性能源供给下，从任意初始条件出发都能产生稳定的周期性运动。

2.采用庞加莱映像，演示强迫VDP 方程在不同参数下所存在四种吸引子，即周期1吸引子、周期2吸引子、不变环面吸引子和奇怪吸引子。

3.对于强迫VDP 方程，在v 和w 为定值条件下，逐渐增大μ值，将出现周期倍分岔和混沌现象。

主要参考资料：
[1] Steven E. Koonin, 秦克诚译. 计算物理学. 北京：高等教育出版社，1993. [2] 马文淦等. 计算物理学. 合肥：中国科学技术大学出版社，1992. [3] 张志涌. 精通MA TLAB6.5. 北京：北京航空航天大学出版社，2003.
完成期限 指导教师
专业负责人
2012年 2 月 23 日
第1章 概述
1.1 自激震荡
自激振动是一种对科学技术非常有意义而在自然界又广泛存在的非线性振动. 我们知道,对线性阻尼振动系统,严格的周期运动只能由受周期性驱动力作用的受迫振动产生;而对非线性系统,有一种自激振动系统,在非振动的能源供给下,它能产生严格的周期运动,这是人们十分感兴趣的现象。

自激系统是一个非线性的有阻尼的振动系统,在振动过程中伴随有能量损耗,但系统存在一种机制,使能量能够由非振动能源通过系统本身的反馈调节,及时适量地得到补充,从而产生一个稳定的不衰减的周期运动,这样的振动称为自激振动[1]。

自激振动现象是一种普遍现象。

如钟摆、弦乐器以及人的心脏的周期性跳动。

活塞发动机的周期性运动等都是利用这种现象来建立不衰减的周期运动；但有些自激振动是十分有害的，这些现象应该设法避免。

1.2 范·德·波耳（Van der pol ）方程
在此课程设计中，我们主要研究范•德•波耳(V an der pol)方程：
(
)02
02
2022
=+--x dt
dx
x
x dt
x d ωμ
（1-1）
所描述的非线性有阻尼的自激振动系统，其中μ是一个小的正的参量，0x 是常数。

下面简称范•德•波耳方程为VDP 方程。

在VDP 方程中，增加外驱动力t V ωcos 项所得到的方程：
(
)0cos 2
02
2022
=++--t V x dt
dx
x
x dt
x d ωωμ
（1-2）
称强迫VDP 方程，其中外驱动力的振幅、角频率分别是V 和ω。

VDP 方程所描述的系统在非线性能源供给下，从任意初始条件出发都能产生稳定的周期性运动。

而对于强迫VDP 方程，在V 和为 ω定值条件下，μ逐渐增大，将出现周期倍分岔和混浊现象。

第2章 自激振荡系统的周期性运动
2.1 课题解析
自激系统是一个非线性有阻尼的振动系统，在运动过程中伴随有能量损耗，但系统存在一种机制，使能量能够由非振动的能源通过系统本身的反馈 调节，及时适量地得到补充，从而产生一个稳定的不衰减的周期运动，这样的振动称为自激振动。

对VDP 方程，可从机械振动角度理解，()220x x --μ是阻尼系数，它是变化的，如果0x x >， 则阻尼系数为正，系统将受阻尼，能量将逐渐减少，但如果
x x <， 则发生负阻尼，意味着不仅不消耗系统的能量，反而给系统提供能量。

此系统能通过自动的反馈调节，使得在一个振动过程中，补充的能量正好等于消耗的能量，从而系统作稳定的周期振动。

2.2 MATLAB 程序设计及演示
取方程(1-1)中的10=x ，10=ω， 3.0=μ，0.66,，0.85，1.08这四个初始条件进行编程，程序详见附录-程序（1）。

2.3 MATLAB演示结果和分析
根据图3.1显示的结果,我们做出以下结论:给出任一初始条件，通过计算机数值求解可以证明它的相轨道都将趋向于一条闭合曲线，这一条闭合曲线，成为极限环，极限环以外的相轨道向里盘旋，而极限环以内的相轨道则向外盘旋，都趋向极限环，说明不论初始情况如何，系统最终都到达以极限环描述的周期性运动[2]。

图3.1 四种初始条件下的轨道
第3章 强迫VDP 方程在不同参数下的四种吸引子
下面研究强迫VDP 方程的行为，我们同时采用时间历程图、相图、庞加莱映像图来研究系统在不同参数条件下的动力学行为，可以看到存在不同的吸引子，即周期1吸引子，周期2吸引子，不变环面吸引子和奇怪吸引子。

3.1庞加莱映射
为了更清楚地了解运动的形态，庞加莱对连续运动的轨迹用一个截面（叫庞加莱截面）将其横截，那么根据轨迹在截面上穿过的情况，就可以简洁地判断运动的形态，由此所得图像叫庞加莱映像。

在截面图上，轨迹下一次 穿过截面的点1+n x 可以看成前一次穿过的点n x 的一种映射。

)(1n n x f x =+（n=0,1,2,…） （3-1）
这个映射就叫庞加莱映射。

它把一个连续的运动化为简洁的离散映射来研究。

在庞加莱映射中的不动点反映了相空间的周期运动，如果运动是二倍周期的，则庞加莱映射是两个不动点，四倍周期则有四个不动点等。

绘制庞加莱映射是在普通的相平面上进行，它不是像画相轨道那样随时间变化连续地画出相点，而是每隔一个外激励周期（ωπ/2=T ）取一个点，例如取样的时刻可以是t=0,T,2T …相应的相点记为),(000y x P 、),(111y x P 、),(222y x P … 这些离散相点就构成了庞加莱映射[3]。

设x y =1，dt
dx y =
2， 则（1-2）式可化为
21y dt
dy =
t V y y y x dt
dy ωωμcos )(12
0221202---=
（3-2）
3.2 强迫VDP 方程在不同参数下的四种吸引子
对于（3-2）式，
t V y y y x dt
dy ωωμcos )(12
022
12
02---=
取120=x ，120=ω，进行以下数值计算研究：
⑴在85.0=μ，1=V ，44.0=ω条件下，存在周期1吸引子，它的周期等于外激励的周期，代表主谐波运动，如图3.1所示：
图 3.1 周期1吸引子
⑵在02.1=μ，1=V ，44.0=ω条件下，存在周期2吸引子，它的周期等于外激励的整数倍，代表次谐波运动，如图3.2所示：
图 3.2 周期2吸引子
⑶在66.0=μ，1=V ，44.0=ω条件下，存在不变环面吸引子，它代表准周期（拟周期）运动，如图3.3所示
图3.3 不变环面吸引子
⑷08.1=μ，1=V ，44.0=ω条件下存在奇怪吸引子，它代表混浊运动。

如图3.4所示
图3.4 奇怪吸引子
⑸保持V和ω为定值，逐渐增大μ，将显示系统状态演化过程全貌的图。

而前四种情况中，看到的只是μ取4个值的片断情况，图形显示，当μ由0.9连续变化到1.2时，系统运动状态逐渐由周期1过渡到周期2（发生了周期倍分岔）再过渡到混浊状态。

如图3.5所示。

图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:4
[t,y]=ode23('zjzdfun',[0:T/1000:50*T],[4,4],[],u(j),x0,w0,v,w); figure
plot(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;
str{1}='庞加莱截面—周期1吸引子';
str{2}='庞加莱截面—周期2吸引子';
str{3}='庞加莱截面—不变环面吸引子';
str{4}='庞加莱截面—奇怪吸引子';
for j=1:4
[t,y]=ode23('zjzdfun',[0:T/1000:50*T],[4,4],[],u(j),x0,w0,v,w); figure
subplot(2,1,1)
plot(t,y(:,1));
title('位移曲线');
xlabel('x');ylabel('v');
subplot(2,2,3)
plot(y(3000:end,1),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(str{j});
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(1)^2)*y(2)-y(1)*w0^2-v*cos(w*t)];
参考文献
[1]丁文静.自激振动，清华大学出版社，2009，(34)：325-328.
[2] 胡静,彭芳麟,管靖,卢圣治.理论力学中非线性问题的MA TLAB数值解,大学物理COLL EGE PHYSICS 第20卷第10期, 2001
[3]张志涌.精通MA TLAB6.5 北京：北京航空航天大学出版社，2003.。