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

Size: px

Start display at page:

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

芷竺凤
7 years ago
Views:

1 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, ),

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

3 (θ = 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,

4 1. G G 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

5 θ, p, ϕ ϕ 2nπ 2(n + 1)π ϕ, θ, p 3.1 f = 0.8,, θ = 2, dθ/dt = 2 θ = 0.1, dθ/dt = 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.

6 3.3 f = 1.65, f = 1.082, f = 1.088, 4, θ = 0.8, dθ/dt = 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

7 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([ ]) 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

8 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 ;

9 end 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;

10 第五节 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 x x2 = E 2 dt 2 2 所以势能是 µ 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.

11 : (0,0); (1,0) (-1,0) < E < 0, x = ±1 E = 0.1, E = 0,, (0,0) 3. E > 0, x = 0, ±1 E = >>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

5-1.5 1.5]) hold on for j=1:length(d0) d=d0(j); [t,u]=ode45(@dbdfun,.

12 阻尼消耗能量,外部驱动补充能量,系统的运动状态解有周期解或混沌解为了掌握运动的整体情况,先画系统的终态解随阻尼系数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([ ]) 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 周期解的情形在map图上周期1吸引子是一个点在频谱图是一个频率 k = 混沌解在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.$ $001]); figure plot(t,u(:,1), r,t1,u1(:,1), g ) xlabel( ); ylabel( ); title( ); %\fs{ d=1.$

13 map k = 1.15 v 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,v ]); figure plot(t,u(:,1), r,t1,u1(:,1), g ) xlabel( ); ylabel( ); title( ); %\fs{ d=1.5, 1 ; d= ; 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);

14 figure set(gcf, unit, normalized, Position,[ ]); 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([ ]) 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([ ]) title( ); xlabel( /Hz ); ylabel( /w ); subplot(2,3,5) %\fs{ } plot(u(2000:300:30000,1),u(2000:300:30000,2), r. ); axis([ ]) title(str{j});

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

. () ; () ; (3) ; (4).. () : P.4 3.4; P. A (3). () : P. A (5)(6); B. (3) : P.33 A (9),. (4) : P. B 5, 7(). (5) : P.8 3.3; P ; P.89 A 7. (6) : P.

() * 3 6 6 3 9 4 3 5 8 6 : 3. () ; () ; (3) (); (4) ; ; (5) ; ; (6) ; (7) (); (8) (, ); (9) ; () ; * Email: [email protected] . () ; () ; (3) ; (4).. () : P.4 3.4; P. A (3). () : P. A (5)(6); B. (3) :