, 10, (Poincare) dθ, ( ) 2 1 dθ cos θ = E 2 dt K V V = cos θ E

1. l, m, θ,, r, ω d F ml d2 θ dt 2 + rldθ dt + mg sin θ = F cos ω dt d 2 θ dt + r dθ 2 m dt + g l sin θ = F ml cos ω dt ω 0 = g/l, mω0, 2 mlω0 2 MLT β = r/2mω 0, f = F/mlω 2 0 = F/mg, ω 0 /ω 0 = 1, ω = ω d /ω 0 d 2 θ dt + 2βdθ + sin θ = f cos ωt 2 dt,β, f, ω, t [ ] ω0 = T 1 ω 0 T 1, [ ] ωd = T 1 ω 0 T 1, [ r mω 0 ] = MT 1 MT 1, [ ] F = mg MLT 2 MLT 2 1. β = 0, f = 0, d2 θ dt 2 + sin θ = 0 1.1 ( 0, ),

, 10, 10 180 180 1.2 (Poincare) dθ, ( ) 2 1 dθ cos θ = E 2 dt K V V = cos θ E

(θ = 0, dθ/dt = 0), K = V = E0, ( ) ( ) (θ = ±π, dθ/dt = 0) θ = π, dθ/dt = 0 θ = π, dθ/dt = 0 E < max(v ), E max(v ),,, 1.3 θ 2π ±π G G 2. β 0, f = 0, d 2 θ dt + 2βdθ 2 dt + sin θ = 0,

1. G G 2. 3. G G ( 2π ) (+2π) 3. β 0, f 0( β = 1/4, ω = 2/3), d 2 θ dt + 2βdθ + sin θ = f cos ωt 2 dt dθ dt = p dp dt = 2β dθ dt dϕ dt = ω sin θ + f cos ωt

θ, p, ϕ ϕ 2nπ 2(n + 1)π ϕ, θ, p 3.1 f = 0.8,, θ = 2, dθ/dt = 2 θ = 0.1, dθ/dt = 0.2 3.2 f = 1.03, θ = 0.1, dθ/dt = 2 θ = 0.8, dθ/dt = 2

3.3 f = 1.65, f = 1.082, f = 1.088, 4, θ = 0.8, dθ/dt = 2 3.4 function dbyd global a f u u=2/3; a=0.5; ZQ=3*pi; f=1.089; [T, Y]=ode45(@dby,[0:ZQ/200:500*ZQ],[-0.8,2]); figure plot(y(31000:end,1),y(31000:end,2)) figure

for j=40001:200:length(y(1:end,1))-1 xx=[xx,y(j,1)]; yy=[yy,y(j,2)]; end plot(xx,yy,.r ) function ydot=dby(t,y) global a f u ydot=[y(2); -sin(y(1))- a*y(2) + f*cos(u*t) ]; 3.5 function djddb figure axis([-8 8-2 2]) hold on % plot([4.5,5.2],[0.8,0.8], g,[4.5,5.2],[0,0],... r,[4.5,5.2],[-0.8~~-0.8], b ); text(5.3,0.8, E<2mgl ); text(5.3,0, E=2mgl ); text(5.3,-0.8, E>2mgl ); xlabel( ); ylabel( d /dt ); % ydot=inline( sqrt(abs(e-1+cos(x))), x, E ); e=[3, 2.5, 2, 1.5,1, 0.5, 0.3, 0.1]; % for k=1:8 if k>3 % E<2mgl

Q{k}=acos(1-e(k)); X=linspace(-Q{k},Q{k},300); y=ydot(x,e(k)); plot(x,y, g,x,-y, g ) elseif k==3 % E=2mgl X=linspace(-2*pi,2*pi,300); y=ydot(x,e(k)); plot(x,y, r,x,-y, r ) else % e>2mgl X=linspace(-2*pi,2*pi,300); y=ydot(x,e(k)); plot(x,y, b,x,-y, b ) end end hold off % [t1,w1]=ode45(@f,[0:0.001:6],[pi/7,0],[]); [t2,w2]=ode45(@f,[0:0.001:6],[pi/3,0],[]); % figure plot(t1,w1(:,1),t2,w2(:,1)); xlabel( ); ylabel( ); legend(, ); % theta=linspace(pi/360,pi-0.1,40); T=[]; options=odeset( Events,@events);% % for i=1:40 ;

end [t,u]=ode45(@f,[0:0.001:20],[theta(i),0],options); T=[T,2*t(end)]; figure plot(theta,t) title( ); xlabel( ); ylabel( ); %----------------------------- function ydot=f(t,y) ydot=[y(2); -9.8*sin(y(1))]; function [value,isterminal,direction]=events(t,y) value=y(2); isterminal=1; direction=1;

第五节 1. 倒摆实验 倒摆与杜芬方程 1.1 倒摆实验演示 1.2 倒摆的简化模型与运动方程 倒摆可以简化成右图的模型,它的运动 可以用杜芬方程描述 d2x dx 3 x + x = f cos ωt + k dt2 dt 改变运动阻尼 可以演示运动状态从 周期解到混沌的变化 3. 杜芬(Duffing)方程 下面用波形图 相图 频谱图和庞加莱截面图 map图 研究系统的运 动 d2x 3.1 无阻尼无驱动情形 x + x3 = 0 2 dt 积分得 µ 2 µ 1 dx 1 1 4 + x x2 = E 2 dt 2 2 所以势能是 µ 1 1 4 V = x x2 2 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 3.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 ) d 2 x dt 2 + kdx dt x + x3 = 0 : (x = ±1 ) 3.2 d 2 x dt 2 + kdx 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.5 1.5-1.5 1.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 3.3.2 混沌解 在map图上 周期2吸引子是两个 点 在频谱图是两个频率 k = 1.35

map k = 1.15 v 0 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.04 0.04 0.94 0.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([-2 2-1.5 1.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([0 0.6 0 0.15]) title( ); xlabel( /Hz ); ylabel( /w ); subplot(2,3,5) %\fs{ } plot(u(2000:300:30000,1),u(2000:300:30000,2), r. ); axis([-2 2-1.5 1.5]) title(str{j});

subplot(2,3,6)%\fs{ } h=plot([0,sin(x0)],[0,cos(x0)], o-, erasemode, xor ); axis([-1 1-1 1]) 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)];