ppt_10_倒摆与杜芬方程

第十讲
倒摆与杜芬方程
1.倒摆实验
1.1倒摆实验演示
1.2倒摆的简化模型与运动方程
倒摆可以简化成右图的模型,它的运动可以用杜芬方程描述，
d 2x dt 2
+k dx dt
−
x +x 3
=f cos ωt 改变运动阻尼，可以演示运动状态从周期解到混沌的变化。

3.杜芬(Duffing)方程
下面用波形图，相图，频谱图和庞加莱截面图（map 图）研究系统的运动。

3.1无阻尼无驱动情形d 2x dt 2−x +x 3
=0积分得
12 dx dt 2+12 12x 4
−x 2 =E
所以势能是
V =
12 12
x 4
−x 2 这时有三个平衡点，x =0(不稳定平衡点),x =±1(是稳定平衡点)。

相平面上奇点:鞍点(0,0);中心型奇点(1,0)和(-1,0)。

根据能量分成三种运动。

1.−0.25<E<0,绕x=±1两个闭轨道，代表在稳定平衡位置附近的周期振动动。

（图中E=−0.1,−0.2）
2.E=0,对应∞型轨道，是同宿轨道,同宿点(0,0)。

3.E>0,绕x=0,±1三个平衡位置的闭轨道，也是周期运动。

图中E=0.2。

>>ezplot(’y^2+x^4/2-x^2’,[-1.6,1.6])
>>hold on
>>ezplot(’y^2+x^4/2-x^2+0.2’)
>>ezplot(’y^2+x^4/2-x^2+0.4’)
>>ezplot(’y^2+x^4/2-x^2-0.4’)
3.2有阻尼无驱动情形
这时的方程成为
d2x dt2+k
dx
dt
−x+x3=0
奇点:原点是鞍点，坐标(x=±1)处两不动点，是吸引子。

整个相平面被分隔成两个区域，不同区的相点分别流向这两个不动点。

3.2有阻尼有驱动情形
这时的方程成为
d2x dt2+k
dx
dt
−x+x3=f cosωt
阻
尼消耗能量,外部驱动补充能量,系统的运动状态解有周期解或混沌解。

为了掌握运动的整体情况,先画系统的终态解随阻尼系数k 变化的分岔图如下。

（0.5≤k ≤1.5,f =1,ω=1.）
function dbd global d
x0=0.1;v0=0.1;d0=0.5:0.002:1.5;
axis([0.51.5-1.51.5])hold on
for j=1:length(d0)d=d0(j);
[t,u]=ode45(@dbdfun,...
[0:2*pi/60:60*pi],[x0,v0]);plot(d,u(901:60:1800,2),’r.’);end
3.3.1周期解的情形
在map 图上，周期1吸引子是一个点，在频谱图是一个频率。

k =1.5
在map 图上，周期2吸引子是两个点，在频谱图是两个频率。

k =1.35
3.3.2混沌解
在map图上的奇怪吸引子，它的频
谱是连续的。

k=1.15
混沌状态下v0相差0.001时形成的
两条分开的曲线。

程序如下
function db global d
%\fs{设v0有微小的变化,比较解的变化情况}
x0=0.1;v0=0.1;d=0.78;
[t,u]=ode45(@dbfun,[0:0.01:100],[x0,v0]);
[t1,u1]=ode45(@dbfun,[0:0.01:100],[x0,v0-0.001]);
figure
plot(t,u(:,1),’r’,t1,u1(:,1),’g’)
xlabel(’时间’);
ylabel(’摆角’);
title(’混沌状态下初条件有微小差异会形成的两条分开的曲线’);
%\fs{当d=1.5,为周期1吸引子;当d=1.35为周期2吸引子;当d=1.15为奇怪吸引子,}
%\fs{读者可以改变d值,以观察不同的情况}
d0=[1.5,1.35,1.15];
str{1}=’庞加莱截面―周期1吸引子’;
str{2}=’庞加莱截面―周期2吸引子’;
str{3}=’庞加莱截面―奇怪吸引子’;
for j=1:3
d=d0(j);
[t,u]=ode45(@dbfun,[0:2*pi/300:200*pi],[x0,v0]);
figure
set(gcf,’unit’,’normalized’,’Position’,[0.040.040.940.8]);
subplot(2,2,1)%\fs{位移曲线}
plot(t,u(:,1))
title(’位移曲线’);
axis([0,150,-2.5,2.5]);
xlabel(’x’);ylabel(’t’);
subplot(2,2,2)%\fs{相图(奇怪吸引子)}
plot(u(20000:end,1),u(20000:end,2))
title(’相图’);
axis([-22-1.51.5])
xlabel(’x’);ylabel(’v’);
Y=fft(u(:,1));%\fs{傅里叶功率分析}
Y(1)=[];n=length(Y);m=fix(n/2);
power=abs(Y(1:m)).^2/n^2;%\fs{功率}
freq=100*(1:n/2)./n;%\fs{频率}
subplot(2,3,4)
plot(freq,power)
axis([00.600.15])
title(’功率谱’);
xlabel(’频率/Hz’);ylabel(’功率/w’);
subplot(2,3,5)%\fs{庞加莱截面}
plot(u(2000:300:30000,1),u(2000:300:30000,2),’r.’);
axis([-22-1.51.5])
title(str{j});
subplot(2,3,6)%\fs{实物模拟图}
h=plot([0,sin(x0)],[0,cos(x0)],’o-’,’erasemode’,’xor’); axis([-11-11])
title(’倒摆运动模拟’);
for i=25000:30000
set(h,’xData’,[0,sin(u(i,1))],’yData’,[0,cos(u(i,1))]);
drawnow
end end
%---------------------
function ydot=dbfun(t,y)
global d
r=1;w=1;
ydot=[y(2);
-y(1)^3+y(1)-d*y(2)+r*cos(w*t)];。